Orthogonal Series Density Estimation for Complex Surveys
Shangyuan Ye, Ye Liang and Ibrahim A. Ahmad

arXiv:1709.06588v2 [stat.ME] 26 Sep 2017

Department of Statistics, Oklahoma State University

Abstract
We propose an orthogonal series density estimator for complex surveys, where samples are neither independent nor identically distributed. The proposed estimator is
proved to be design-unbiased and asymptotically design-consistent. The asymptotic
normality is proved under both design and combined spaces. Two data driven estimators are proposed based on the proposed oracle estimator. We show the efficiency of
the proposed estimators in simulation studies. A real survey data example is provided
for an illustration.
Keywords: Nonparametric, asymptotic, survey sampling, orthogonal basis, HorvitzThompson estimator, mean integrated squared error.

1

Introduction

Nonparametric methods are popular for density estimations. Most work in the area of
nonparametric density estimation was for independent and identically distributed samples.
However, both assumptions are violated if the samples are from a finite population using
a complex sampling design. Bellhouse and Stafford (1999) and Buskirk (1999) proposed
kernel density estimators (KDE) by incorporating sampling weights, and their asymptotic
properties were studied by Buskirk and Lohr (2005). Kernel methods for clustered samples
and stratified samples were studied in Breunig (2001) and Breunig (2008), respectively.
One disadvantage of the KDE is that all samples are needed to evaluate the estimator.
However, in some circumstances, there is a practical need to evaluate the estimator without
using all samples for confidentiality or storage reasons. For example, many surveys are
routinely conducted and sampling data are constantly collected. Data managers want to
publish exact estimators without releasing all original data. In Section 6, we provide a real
data example from Oklahoma M-SISNet, which is a routinely conducted survey on climate
policies and public views. The orthogonal series estimators are useful alternatives to KDEs,
without needing to release or store all samples.
1

The basic idea of the orthogonal series method is that any square integrable function
f , in our case a density function, can be projected onto an orthogonal basis {ϕj }: f (x) =
P∞
j=0 θj ϕj (x), where
θj =

Z

ϕj (x)f (x)dx = E(ϕj (X))

(1)

is called the jth Fourier coefficient. Some of the work in orthogonal series density estimation (OSDE) was covered in monographs by Efromovich (1999) and Tarter and Lock
(1993), among others. Efromovich (2010) gave a brief introduction of this method. Walter
(1994) discussed properties of different bases. Donoho et al. (1996) and Efromovich (1996)
studied data driven estimators. Asymptotic properties were studied by Pinsker (1980) and
Efromovich and Pinsker (1982).
In this paper, we study the OSDE for samples from complex surveys. To the best of
our knowledge, no previous work has been done on developing OSDE for finite populations.
We propose a Horvitz-Thompson type of OSDE, incorporating sampling weights from the
complex survey. We show that the proposed OSDE is design-unbiased and asymptotically
design-consistent. We further prove the asymptotic normality of the proposed estimator. We
compare the lower bound of minimax mean integrated squared error (MISE) with the I.I.D.
case in Efromovich and Pinsker (1982). We propose two data driven estimators and show
their efficiency in a simulation study. Finally, we analyze the M-SISNet survey data using
the proposed estimation. All proofs to theorems and corollaries are given in the appendix.

2

Notations

Consider a finite population labeled as U = {1, 2, ..., N}. A survey variable x is associated
with each unit in the finite population. A subset s of size n is selected from U according to
some fixed-size sampling design P̧(·). The first and second order inclusion probabilities from
the sampling design P̧(·) are πi = Pr(i ∈ s) and πij = Pr(i, j ∈ s), respectively. The inverse
of the first order inclusion probability defines the sampling weight di = πi−1 , ∀i ∈ s.
The inference approach used in this paper for complex surveys is the combined designmodel-based approach originated in Hartley and Sielken (1975). This approach accounts
for two sources of variability. The first one is from the fact that the finite population is
a realization from a superpopulation, that is, the units xU = {x1 , x2 , ..., xN } are considered independent random variables with a common distribution function F , whose density
function is f . The second one is from the complex sampling procedure which leads to a
2

sample x = {x1 , x2 , . . . , xn }. Denote w = {w1 , w2 , . . . , wn } design variables that determine
the sampling weights. The sampling design P̧(·) is embedded within a probability space
(S, J̧, PP̧ ). The expectation and variance operator with respect to the sampling design are
denoted by EP̧ (·) = EP̧ (· | xU ) and VarP̧ (·) = VarP̧ (· | xU ), respectively. The superpopulation ξ, from which the finite population is realized, is embedded within a probability space
(Ω, F̧, Pξ ). The sample x and the design variables w are ξ-measurable. The expectation
and variance operator with respect to the model are denoted by Eξ (·) and Varξ (·), respectively. Assume that, given the design variables w, the product space, which couples the
model and the design spaces, is (Ω × S, F̧ × J̧, Pξ × PP̧ ). The combined expectation and
variance operators are denoted by EC (·) and VarC (·), where EC (·) = Eξ [EP̧ (· | xU )] and
VarC (·) = Eξ [VarP̧ (· | xU )] + Varξ [EP̧ (· | xU )].

3

Main Results

Consider a sample s = {x1 , x2 , ..., xn } drawn from a finite population xU using some fixedsize sampling design P̧(·). Our goal is to estimate the hypothetical density function f of the
superpopulation. Equation (1) implies that θj can be estimated using the Horvitz-Thompson
(HT) estimator for the finite population mean
θ̂j = N

−1

n
X

di ϕj (xi ),

(2)

i=1

where N is the finite population size and di = πi−1 is the sampling weight for unit i. The
HT estimator is a well known design unbiased estimator (Fuller, 2009). The basis {ϕj } can
be Fourier, polynomial, spline, wavelet, or others. Properties of different bases are discussed
in Efromovich (2010). We consider the cosine basis throughout the paper, which is defined
√
as {ϕ0 = 1, ϕj = 2 cos(πjx)}, j = 1, 2, · · · , x ∈ [0, 1]. Regarding the compact support [0, 1]
for the density, we adopt the argument in Wahba (1981):“it might be preferable to assume
the true density has compact support and to scale the data to interior of [0, 1].” Analogous
to Efromovich (1999), we propose an orthogonal series estimator in the form
fˆ(x) = fˆ(x, {wj }) = 1 +

∞
X

wj θ̂j ϕj (x),

(3)

j=1

where θ̂j is the HT estimator for the Fourier coefficient as in (2) and wj ∈ [0, 1] is a shrinking
R1
coefficient. Note that θ0 = 0 f (x)dx = 1. If xU is known for all units in the finite population,
3

we can write the population estimator for f (x) as
fU (x) = fU (x, {wj }) = 1 +
where θU,j = N

PN
−1

i=1

∞
X

wj θU,j ϕj (x),

(4)

j=1

ϕj (xi ).

The following theorems and a corollary show properties of our proposed estimator under
both design and combined spaces. Theorem 1 considers unbiasedness and consistency under
the design space.
PP
P∞ 2
πik
Theorem 1 Suppose f ∈ L2 (R), δ = N −1
i=1 wi < ∞.
i6=k πi πk − N < ∞ and
Then, the estimator fˆ(x, {wj }) is design-unbiased and asymptotically design-consistent for
fU (x, {wj }), i.e.,
h
h
i
i
ˆ
ˆ
EP̧ f (x, {wj }) = fU (x, {wj }) and ΓP̧ = VarP̧ f (x, {wj }) → 0 as N → ∞.

Theorem 2 shows the asymptotic normality of the proposed estimator fˆ(x, {wj }) under
the design space.
Theorem 2

Suppose that all assumptions in Theorem 1 hold. As N → ∞,
fˆ(x, {wj }) − fU (x, {wj })
Γ̂P̧

L
P̧
−−→ N(0, 1),

(5)

where
Γ̂P̧ = N −1

J
X

wj2 (1 + 2−1/2 θ̂2j + δ θ̂j2 )(1 + 2−1/2 ϕ2j (x)).

j=1

We then show the asymptotic normality of the proposed estimator fˆ(x, {wj }) under the
combined inference. Define a Sobolev Class of k-fold differentiable densities as F̧(k, Q) =
P∞
P
2k 2
{f : f (x) = 1 + ∞
j=1 (πj) θj ≤ Q < ∞}, k ≥ 1. Note that for any
j=1 θj ϕj (x),
f ∈ F̧(k, Q), f is 1-periodic, f (k−1) is absolute differentiable and f (k) ∈ L2 (R).
Theorem 3

Suppose that f ∈ F̧(k, Q) and all assumptions in Theorem 2 hold. Then,
fˆ(x, {wj }) − f (x) LC
h
i −→ N(0, 1) as N → ∞,
VarC fˆ(x, {wj })

(6)

h
i
P
where VarC fˆ(x, {wj }) = N −1 Jj=1 wj2 bj (1 + 2−1/2 ϕ2j (x)) and bj = 2 + 21/2 θ2j + (δ − 1)θj2 +
oN (1).
4

The following corollary is a direct result of using Theorem 3 and Efromovich and Pinsker
(1982). It shows the lower bound of the minimax MISE for the proposed estimator fˆ(x, {wj })
under the Sobolev class.
Corollary 1 Let f ∈ F̧(k, Q) and fˆ(x, {wj }) be the estimator in Theorem 3. The lower
bound of the minimax MISE, under the combined inference approach, is given by:
h
i
R(F̧) = inf sup MISEC fˆ(x, {wj }) ≥ N −2k/(2k+1) P (k, Q, b)(1 + oN (1)),
(7)
{wj }
f ∈F̧(k,Q)
o2k/(2k+1)
n
k
and b = 2.
where P (k, Q, b) = Q1/(2k+1) π(k+1)b

Remark that this lower bound is of the same form as the I.I.D. case in Efromovich and Pinsker
(1982), but with b = 2 instead of b = 1.

4

Data Driven Estimators

The choice of shrinking coefficients ŵj is not unique. To get a proper data driven estimator,
we start with the oracle estimator (3), and then obtain ŵj by minimizing the MISE for the
oracle estimator. Here, we propose two estimators: a truncated estimator and a smoothed
truncated estimator, mimicking those in the I.I.D. case.
The truncated estimator, denoted by fˆT , is an estimator with ŵj = 1 for j ≤ J, and

ŵj = 0 for j > J. Alternatively we can write ŵj = Ij≤J . Then, only the truncation
parameter J needs to be estimated. Notice that the MISE of this estimator is
J h
h
i Z
i X
2
ˆ
MISEC f (x, {wj }) =
VarC (θ̂j ) − θj − f 2 (x)dx.
j=1

Since f 2 (x)dx is fixed and an unbiased estimator for θj2 is θ̂j2 −N −1 bj , a data-driven estimate
for J can be obtained from
R

Jˆ = arg min

J
X
(2N −1 b̂j − θ̂j2 ),
j=1

where b̂j is the plug-in estimator of bj . In practice, the solution is obtained through a
numerical search. Efromovich (1999) suggests to set the upper bound for Jˆ to be ⌊4+0.5 ln n⌋
for the search. Theoretically, the minimum of the MISE can be approximated in the following
corollary.
5

Corollary 2

Let f ∈ F̧(k, Q), k > 1/2. The MISE of fˆT is minimized when
J ≈ N 1/(2k+1) H1 (k, b, c),

(8)

and the minimum is approximately
R(fˆT ) = MISEC (fˆT (x, {ŵj })) ≈ N −2k/(2k+1) H2 (k, b, c),
−1/(2k+1)

where H1 (k, b, c) = b
and c is a constant.



2k+1
(2k+2)c

−1/(2k+1)

2k/(2k+1)

, H2 (k, b, c) = b

(9)


2k+1
(2k+2)c

−1/(2k+1)

,

One possible modification for fˆT is to shrink each Fourier coefficient toward zero. We call
this estimator the smoothed truncated estimator, denoted by fˆS . It is constructed similarly
as the truncated estimator, with the first J Fourier coefficients shrunk by multiplying the
optimal smoothing coefficients wj∗, obtained from the proof of Corollary 1. Mathematically,
ŵj = ŵj∗ Ij≤J , where ŵj∗ = (θ̂j2 − N −1 b̂j )/θ̂j2 is a direct plug-in estimator for wj∗ .

A potential problem of the nonparametric density estimation is that the estimator may
not be a valid density function. A simple modification is to define the L2 -projection of fˆT
(or fˆS ) onto a class of non-negative densities, f˜T (x) = max{0, fˆT (x) − const.}, where the
normalizing constant is to make f˜T integrate to 1. It has been proved that the constant
always exists and is unique (Glad et al., 2003).

5

Simulation

We compared our proposed estimators with the series estimator that ignores the finite population and sampling designs, through a Monte Carlo simulation study. We considered
estimating density functions for three sampling designs: (1) the simple random sample without replacement (SRSWOR), (2) the stratified sampling and (3) the Poisson sampling. Note
that the Poisson sampling has a random size with units independently sampled and hence
violates our assumption of fixed size sampling.
1. For the SRSWOR, we considered two superpopulations: the standard normal distribution N(0, 1) and a mixture normal distribution 0.4N(−1, 0.5) + 0.6N(1, 1).
2. For the stratified sampling, we considered two superpopulations: a two-component
mixture normal 0.4N(−1, 0.5) + 0.6N(1, 1) and a three-component mixture normal
6

0.3N(−1, 0.15) + 0.4N(0, 0.15) + 0.3N(1, 0.15). We designed two strata for the twocomponent mixture and three strata for the three-component mixture. A proportional
stratified sampling is used.
3. For the Poisson sampling, we considered the same two superpopulations as in (1).
We specified the expected sample size for the Poisson sampling to be n, and generated the first order inclusion probabilities for the Poisson sampling using the function
“inclusionprobabilities” in the R package “sampling” (Till and Matei, 2016).
For all cases, we considered a finite population of size N = 1, 000 drawn from each of the
superpopulations. We repeated drawing the finite population for m1 = 100 times. For each
of the finite population, we drew samples according to the sampling design, with increasing
sample sizes: n = 20, 40, 60 and 80. The replication number for each finite population is
m2 = 10, 000. The performance of estimators is measured by a Monte Carlo approximation
of the MISE:
Z
m1 X
m2 h
i2
1 X
f˜ij (x) − f (x) dx.
MISEMC (f˜) =
mm
1

2 i=1 j=1

The results of the simulation study are shown in Table 1. In general, the I.I.D. series
estimator, which ignores the sampling design, performs the worst in all cases. However,
it is not surprising to see that the improvement for the proposed estimators is much more
in stratified sampling than in SRSWOR or Poisson sampling. It confirms the necessity of
incorporating stratification sampling weights into the series estimator for a complex survey.
Lastly, the smoothed truncated estimator performs better than the truncated estimator in
most scenarios.

6

Oklahoma M-SISNet Survey

The Oklahoma Weather, Society and Government Survey conducted by Meso-Scale Integrated Sociogeographic Network (M-SISNet) measures Oklahomans’ perceptions of weather
in the state, their views on government policies and societal issues and their use of water and
energy. The survey is routinely conduced at the end of each season. Until the end of 2016, 12
waves of survey data have been collected. It is desired that estimates can be obtained without constantly pulling out the original data. The sampling design has two separated phases.
In Phase I, a simple random sample of size n = 1, 500 is selected from statewide households.
In Phase II, a stratified oversample is selected from five special study areas: Payne County,
7

Oklahoma City County, Kiamichi County, Washita County and Canadian County. In each
stratum, the sample size is fixed to be 200. The second phase can be viewed as a stratified
sampling over the entire state with six strata: n1 = · · · = n5 = 200 and n6 = 0, where the
sixth stratum contains households not in the five special study areas. This design with oversampling is not a typical fixed-size complex survey. The first-order inclusion probabilities
are approximately πhi = nh /Nh + n/N, for i = 1, . . . , Nh and h = 1, . . . , 6. Note that for
units not in the five areas, this inclusion probability is simply n/N. We presents OSDEs for
two continuous variables for illustration: the monthly electricity bill and the monthly water
bill. Figure 1 shows OSDEs of the two variables for all seasons in 2015.

Acknowledgement
This research is partially supported by National Science Foundation under Grant No. OIA1301789.

Appendix
Proof of Theorem 1
Proof. We first show that fˆ(x, {wj }) is design-unbiased:


∞
h
i
X
EP̧ fˆ(x, {wj }) = EP̧ 1 +
wj θ̂j ϕj (x)
j=1

= 1+
= 1+

∞
X

j=1
∞
X

wj EP̧ (θ̂j )ϕj (x)
wj θU,j ϕj (x)

j=1

= fU (x, {wj }).
It remains to show that fˆ(x, {wj }) is asymptotically design-consistent, that is, the design-variance
of fˆ(x, {wj }) approaches zero in the limit. We need the simple fact that
√
ϕ2j (x) = [ 2 cos(πjx)]2 = 1 + cos(π2jx) = 1 + 2−1/2 ϕ2j (x).

8

Then, we have


ΓP̧ = VarP̧ 1 +
=

∞
X

∞
X
j=1



wj θ̂j ϕj (x)

wj2 ϕ2j (x)VarP̧ (θ̂j )

j=1

=

∞
X

wj2

j=1

h

#
" n
i
X
di ϕj (xi ) ,
1 + 2−1/2 ϕ2j (x) N −2 VarP̧
i=1

and
VarP̧

"

n
X
i=1

#

di ϕj (xi )

= VarP̧
= EP̧
=

N
X
i=1

−
=

"

"

N
X

i=1

Ii di ϕj (xi )

i=1

N
X

#2

Ii di ϕj (xi )

i=1

−

(

EP̧



EP̧ (Ii2 )d2i EP̧ ϕ2j (xi ) +

(N
X

N
X

#

EP̧ (Ii )di EP̧ [ϕj (xi )]

i=1

"

N
X

#)2

Ii di ϕj (xi )

i=1

XX

πik di dk EP̧ [ϕj (xi )] EP̧ [ϕk (xk )]

i6=k
)2

h
i XX π
ik 2
2
EP̧ 1 + 2−1/2 ϕ2j (xi ) +
θ − N 2 θU,j
πi πk U,j
i6=k

2
= N (1 + 2−1/2 θU,2j + δθU,j
)

≤ N M,
2 ≤ M < ∞ for every j.
where 1 + 2−1/2 θU,2j + δθU,j
P
2
Hence, ΓP̧ ≤ N −1 M ∞
j=1 wj → 0 as N → ∞.

Proof of Theorem 2
Proof. By the definition of θ̂j and θU,j , we have
fˆ(x, {wj }) = 1 +
= 1+

∞
X

wj θ̂j ϕj (x)

j=1

N
X
i=1

9

Ii di

∞
X
j=1

wj ϕj (x)ϕj (xi ),

and


E Ii di

∞
X
j=1



wj ϕj (x)ϕj (xi ) =
=

∞
X

wj ϕj (x)E [ϕj (xi )]

j=1

∞
X

wj θU,j ϕj (x).

j=1

Also, from the proof of Theorem 1, we have


∞
∞
X
X
2
wj2 (1 + 2−1/2 θU,2j + δθU,j
)
wj ϕj (x)ϕj (xi ) =
Var Ii di
j=1

j=1

≤ B

∞
X
j=1

wj2 < ∞ by assumption.

Therefore, by the Lindeberg-Lévy central limit theorem, we have
fˆ(x, {wj }) − fU (x, {wj }) LP̧
−−→ N (0, 1).
ΓP̧

(10)

It remains to show that Γ̂P̧ is consistent for ΓP̧ under design, or equivalently,
P
P̧
|Γ̂P̧ − ΓP̧ | −−→ 0, as n → N.

(11)

Condition (11) can be proved by using the facts that θ̂j is design unbiased and E(θ̂j2 ) = θj2 +
Var(θ̂j ) → θj2 as n → N .

Then, Theorem 2 is proved by using the equations (10) and (11) in conjunction with Slutsky’s
theorem.

Proof of Theorem 3
Proof. Since fU (x, {wj }) is the standard OSDE from an I.I.D. sample which is the finite population,
then
fU (x, {wj }) − f (x) Lξ
−→ N (0, 1).
Varξ [fU (x, {wj })]

10

(12)

The asymptotic distribution of the I.I.D. OSDE under Sobolev class is obtained from Efromovich
(1999), Chapter 7. Also,
J
h
i
h
i
X
VarC wj θ̂j ϕj (x)
VarC fˆ(x, {wj }) =
j=1

=

J
X

wj2 (1 + 2−1/2 ϕ2j (x))VarC (θˆj )

(13)

j=1

Next, we calculate the variance of θ̂j by using Theorem 1:
h
i
h
i
VarC (θ̂j ) = Eξ VarP̧ (θ̂j ) + Varξ EP̧ (θ̂j )
h
i
2
) + Varξ (θU,j )
= Eξ N −1 (1 + 2−1/2 θU,2j + δθU,j
h
i
2
) + Varξ (θU,j )
= N −1 1 + 2−1/2 θ2j + δEξ (θU,j

(14)

2 ) and Var (θ
Then, we evaluate Eξ (θU,j
ξ U,j ) separately. Based on a standard result in the I.I.D. case,
we have

Varξ (θU,j ) = N −1 (1 + 2−1/2 θ2j − θj2 )

(15)

and
2
Eξ (θU,j
) = E2ξ (θU,j ) + Varξ (θU,j )

= N −1 (1 + 2−1/2 θ2j − θj2 ) + θj2 .
Then, plug equations (15) and (16) into (14), we have
h
i
VarC (θ̂j ) = N −1 2 + 21/2 θ2j + (δ − 1)θj2 + oN (1) = N −1 bj .

(16)

(17)

Hence, plug (17) into (13) we can get the variance of fˆ under the combined inference approach.
Finally, apply Theorem 5.1 in Bleuer and Kratina (1999), Theorem 3 is proved.

Proof of Corollary 1
Proof. The proof is similar to Efromovich and Pinsker (1982). We sketch the steps as follows. We
first evaluate the linear minimax MISE for the functions in the Sobolev class defined above. That
is, we optimize wj∗ ’s that minimize MISEC (fˆ). Notice that EC (θ̂j ) = Eξ [EP̧ (θ̂j )] = Eξ (θU,j ) = θj

11

implying that θ̂j is an unbiased estimator of θj . Therefore,
Z

h
i
MISEC fˆ(x, {wj }) = EC
(f − fˆ)2

∞ n
o
i
h
X
wj2 VarC (θ̂j ) + θj2 − 2wj θj2 + θj2 .
=

(18)

j=1

A straightforward calculation yields that
wj∗ =

θj2
θj2 + VarC (θ̂j )

.

(19)

Plug equation (19) into (18),
h
i
sup MISEC fˆ(x, {wj })
{wj } f ∈F̧(k,Q)
∞
X
θj2 VarC (θ̂j )
≥
sup
,
2
f ∈F̧(k,Q) j=1 θj + VarC (θ̂j )

RL (F̧) =

inf

(20)

where VarC (θ̂j ) is of the form (17). Plug (17) into (20), and use the Lagrange multiplier to show
that the maximum of (6) is attained at
θj2 = N −1 (µ/(πj)k − bj )+ ,
where µ is determined by the constraint
obtain

P∞

2k 2
j=1 (πj) θj

(21)

≤ Q. Plug equation (21) back to (20), we

RL (F̧) ≥ N −2k/(2k+1) P (k, Q, b).
Pinsker (1980) shows that for Sobolev ball F̧, the linear minimax risk is asymptotically equal to
the minimax risk, that is, R(F̧) = RL (F̧)(1 + oN (1)). Therefore Corollary 1 is proved.

Proof of Corollary 2
Proof. Let ŵj = Ij≤J . Plug equation (17) into (18), we have
R(fˆT ) = N −1

J
X
j=1

bj +

∞
X

j=J+1

θj2 ≈ N −1 bJ +

12

∞
X

j=J+1

θj2 .

(22)

Notice that for f ∈ F̧(k, Q). By a straightforward calculation, we have θj2 = cj −2(k+1) (Efromovich,
1999). Therefore,
∞
X

j=J+1

θj2

≈c

Z

∞

j −2(k+1) dj =

J

c
J −2k−1 .
2k + 1

(23)

Plug (23) into (22) and optimize J, Corollary 2 is proved.

References
Bellhouse, D. and Stafford, J. (1999), ‘Density estimation from complex surveys’, Statistica
Sinica 9, 407–424.
Bleuer, S. and Kratina, I. (1999), ‘On the two-phase framework for joint model and designbased inference’, The Annals of Statistics 33, 2789–2810.
Breunig, R. (2001), ‘Density estimation for clustered data’, Econometric Reviews 20, 353–
367.
Breunig, R. (2008), ‘Nonparametric density estimation for stratified samples’, Statistics and
Probability Letters 78, 2194–2200.
Buskirk, T. (1999), Using nonparametric methods for density estimation with complex survey
data, Technical report, PhD thesis, Department of Mathematics, Arizona State University.
Buskirk, T. and Lohr, S. (2005), ‘Asymptotic properties of kernel density estimation with
complex survey data’, Journal of Statistical Planning and Inference 128, 165–190.
Donoho, D., Johnstone, I., Kerkyacharian, G. and Picard, D. (1996), ‘Density estimation by
wavelet thresholding’, Annals of Statistics 24, 508–539.
Efromovich, S. (1996), ‘Adaptive orthogonal series density estimation for small samples’,
Computational Statistics and Data Analysis 22, 599–617.
Efromovich, S. (1999), Nonparametric Curve Estimation: Methods, Theorey and Applications, New York: Springer.
Efromovich, S. (2010), ‘Orthogonal series density estimation’, WIREs Comp Stat 2, 467–476.

13

Efromovich, S. and Pinsker, M. (1982), ‘Estimation of square-integrable probability density
of a random variable’, Problems of Information Transmission 18, 19–38.
Fuller, W. (2009), Sampling Statistics, Wiley, New York.
Glad, I., Hjort, N. and Ushakov, N. (2003), ‘Correction of density estimators that are not
densities’, Scandinavian Journal of Statistics 30, 415–427.
Hartley, H. and Sielken, R. (1975), ‘A super-population viewpoint for finite population sampling’, Biometrics 31, 411–422.
Pinsker, M. (1980), ‘Optimal filtration of square-integrable signals in Gaussian noise’, Problems Inform. Transmission 16, 53–68.
Tarter, M. and Lock, M. (1993), Model-Free Curve Estimation, New York: Chapman and
Hall.
Till, Y. and Matei, A. (2016), sampling: Survey Sampling. R package version 2.8.
URL: https://CRAN.R-project.org/package=sampling
Wahba, G. (1981), ‘Data-based optimal smoothing of orthogonal series density estimates’,
The Annals of Statistics 9, 146–156.
Walter, G. (1994), Wavelets and other Orthogonal Systems with Applications, London: CRC
Press.

14

Table 1: Monte Carlo approximation of MISE for three sampling designs and two superpopulations. The finite population size is N = 1, 000. The replication size of the finite
population is m1 = 100, and the replication size of the sample is m2 = 10, 000. Three
estimators are compared: the truncated estimator, the smoothed estimator and the series
estimator ignoring finite population and sampling design (I.I.D.).
SRSWOR
n
20
40
60
80

Standard Normal
Mixture Normal
Truncated Smoothed I.I.D. Truncated Smoothed
0.0232
0.0220
0.0290
0.0498
0.0480
0.0150
0.0140
0.0157
0.0311
0.0318
0.0116
0.0109
0.0121
0.0226
0.0234
0.0094
0.0089
0.0100
0.0173
0.0180

I.I.D.
0.0535
0.0388
0.0335
0.0219

n
20
40
60
80

Poisson Sampling
Standard Normal
Mixture Normal
Truncated Smoothed I.I.D. Truncated Smoothed
0.0497
0.0481
0.0527
0.0580
0.0442
0.0281
0.0270
0.0392
0.0344
0.0294
0.0241
0.0229
0.0237
0.0283
0.0280
0.0201
0.0190
0.0211
0.0235
0.0234

I.I.D.
0.0705
0.0399
0.0322
0.0285

n
20
40
60
80

Stratified Sampling
Two Strata
Three Strata
Truncated Smoothed I.I.D. Truncated Smoothed
0.0415
0.0409
0.0739
0.2847
0.2826
0.0231
0.0230
0.0688
0.2731
0.2718
0.0181
0.0180
0.0672
0.0426
0.0419
0.0142
0.0142
0.0675
0.0412
0.0406

I.I.D.
0.3106
0.3309
0.1132
0.1175

15

0.020

(a) Water Bill

0.000

0.005

Density
0.010

0.015

Winter
Spring
Summer
Fall

0

50

100
150
200
Amount (in dollars)

250

300

(b) Electricity Bill

0.000

Density
0.004
0.008

Winter
Spring
Summer
Fall

0

50

100
150
200
Amount (in dollars)

250

300

Figure 1: OSDEs of the electricity bill and the water bill for seasonal waves in 2015.

16

