H ∞ Optimal Nonparametric Density Estimation from
Quantized Samples
M. Nagahara1 , K. I. Sato2 , and Y. Yamamoto1
Graduate School of Informatics, Kyoto University,
2
Graduate School of Economics, Kyoto University,
Sakyo-ku Yoshida-Honmachi, Kyoto, 606-8501, Japan
1

Abstract
In this paper, we study nonparametric density estimation from quantized samples. Since quantization
decreases the amount of information, interpolation (orestimation) of the missing information is needed. To
achieve this, we introduce sampled-data H ∞ control
theory to optimize the worst case error between theoriginal probability density function and the estimation. The optimization is formulated by linear matrix
inequalities and equalities. A numerical example is illustrated to show the effectiveness of our method.

1

Introduction

Probability density functions (PDF’s, henceforth) or
probability distributions play a fundamental role inboth analysis and design of stochastic systems. Theoretically, it is often assumed that such information
about randomness is given a priori. However, in practice, we usually have only limited information about
that. On one hand, the random nature of such systems
hinders us to obtain the exact distribution. On theother hand, because the analysis and design are usually executed on a digital processor, acquired data must
be discretized in both time and value and this severely
reduces the information. In this paper, we propose amethod that weakens the latter difficulty. That is, weare concerned about estimating a PDF from quantized
data. Our analysis assumes no specific function form
for the PDF, and hence our method is supposed to beclassified into nonparametric estimation. At the sametime, it might be worth mentioning that we exploit a
couple of ideas which were given rise in the signal processing and control literatures.
The first idea that we rely on owes Widrow [13, 14].
These papers study the estimation problem from quantized samples by applying Shannon’s sampling theorem[8, 10]. We refer to the Widrow’s result as Widrow’squantizetion theorem. This theorem assumes that theoriginal characteristic function is band-limited up to
the frequency π/∆, where ∆ is the step size of theuniform quantizer. This frequency corresponds to the
Nyquist frequency in Shannon’s sampling theorem. The
assumption which is called band-limiting condition isquite restrictive because many random variables suchas uniform, Gaussian, Erlangian, or Rayleigh ones do
not satisfy the assumption [6].
The second helpful idea has come from the control
theory. In order to avoid the stringent assumptionof Shannon’s theorem, sampled-data H ∞ infty control theory proved useful in [5]. The paper shows thatthe sampled-data H ∞ optimization provides a way toobtain a filter without the band-limiting assumption.
By altering the method of [5] into that suitable for
Widrow’s quantization theorem, we will be able to obtain a better estimation.
In this paper, this alternative way to estimating a
PDF is proposed. Unlike Widrow’s method, we areunable to obtain the exact PDF. Such a restriction isinevitable because we dismiss the band-limiting condition. Instead, we apply the sampled-data control theory[2] that enables us to minimize the undesirable error
in the analog domain. A designed filter, which generates an estimated PDF, will be combined with an
up-sampler to compose an interpolator [11]. Like themethod of [5], our estimation can be accomplished by
efficient numerical optimization.
The organization of this paper is as follows. Section
2 defines nonparametric density estimation problem
which we discuss in this paper. In section 3, we introduce Widrow’s quantization theorem, and we point out
problems in Widrow’s estimation. To solve these problems, we propose a new estimation method by sampleddata H ∞ optimization in section 4. Numerical examples, given in section 5, illustrate the effectiveness of our
method compared with Widrow’s estimation. Section
6 concludes our results and shows a future work.
Notations
We use the following notations. R and Z are the setsof the real numbers and the integers, respectively. L2
is the Lebesgue spaces consisting of the square integrable real functions on R. For w ∈ L2 , F w denotes
the Fourier transform of w. We denote by χI the indicator function on a set I ⊂ R; χI (x) = 1 if x ∈ I andχI (x) = 0 if x ∈/ I. For a matrix M , M ⊤ and trM are
respectively the transpose and the trace of M . For asymmetric matrix P , P ≻ 0, P  0 and P ≺ 0 denotepositive, nonnegative, and respectively negative definite matrices. For a linear time-invariant (continuoustime or discrete-time) system with state space matrices

Qx

This problem is ill-posed because there are infinitely
many solutions. To avoid ill-posedness, we have to restrict the function space to which the original PDF belongs. A unique solution to the problem is obtained by
Widrow et al. by restricting the function space to thatof band-limited functions [13, 14]. We introduce thissolution in the next section.

6
4
2
-5

-3

1

-1

3

5

x

0

3

-2
-4
-6

Fig. 1: Uniform quantizer with ∆ = 2 and M = 3

Widrow’s Quantization Theorem

In this section, we introduce Widrow’s quantization
theorem.
First, we analyze the uniform quantization Q (see(1)) of a PDF p(x). The samples drawn from p(x)and quantized by Q can be interpreted as samples ofa discrete PDF c[k]. This discrete PDF is obtained byarea-sampling [13, 14], denoted by Se∆ p, that is,
Z (k+1/2)∆
c[k] =
p(x)dx =: (Se∆ p)[k](2)(k−1/2)∆{A, B, C, D},

A
C

B
D



(s) or



A
C

B
D



(z)

denotes the transfer function.

2

Density Estimation Problem

In this section, we define our problem. Assume that
there are samples drawn from a PDF p(x), which is
unknown. Our concern in the present paper is how
to estimate p(x) from the samples. In ordinary practice the acquired samples are stored and processed indigital systems, and hence they must be converted todigital data, that is quantized. We assume that thisquantization is a uniform one with the step size ∆ andthe no-overload input range ((−M − 1)∆, (M + 1)∆)(see [7]), which is defined by
Q : R → A := {−M ∆, (−M + 1)∆, . . . , M ∆},



k∆,
k − 21 ∆ ≤ x < k + 12 ∆,k = −M + 1, . . . , M − 1, (1)

Qx :=

M ∆,
x ≥ M − 21 ∆,




−M ∆,
x < −M − 21 ∆.

Fig. 1 shows a uniform quantizer with ∆ = 2 andM = 3.
Then, our problem is formulated as follows:
Problem 1 Given samples x
e1 , xe2 , . . . , xeN which are
i.i.d. (independent, identically distributed) and quantized by the uniform quantizer defined by (1), estimate the original PDF p(x) from which the samples are
drawn.

This can be rewritten byZ ∞c[k] =
χ[−∆/2,∆/2) (k∆ − x)p(x)dx−∞= (q∆ ∗ p)(k∆) =: S∆ (q∆ ∗ p)[k],q∆ (x) := χ[−∆/2,∆/2) (x).

Thus, quantization can be interpreted as convolutionby q∆ and ideal sampling S∆ . Note that the Fouriertransform (frequency response) of q∆ is given by2
1  jω 2(3)e ∆ − e−jω ∆ .
Q∆ (jω) =jω
Based on this idea, Widrow et al. proposed a reconstruction scheme by using Shannon’s sampling theorem.
To see this, we define the space BL2 of band-limited
functions:
BL2 := p ∈ L2 : (F p)(ω) = 0, |ω| > π/∆ .

Then, Widrow et al. have shown a solution to Problem 1 as follows [13, 14].
Theorem 1 (Widrow et al.) Assume that p ∈ BL2 .
Then p(x) can be perfectly reconstructed from the discrete PDF of the quantized samples {exn }.
The reconstruction procedure by this theorem is asfollows. First, construct a histogram from the samples {exn }, with its intervals (−∞, −(M + 1/2)∆),[(k − 1/2)∆, (k + 1/2)∆), k = −M + 1, . . . , M − 1, and[(M − 1/2)∆, ∞). If we have sufficiently many samples, {ck } can be approximated from the frequency ofeach interval divided by the number of the samples. Letr = q∆ ∗p. Then, if the band-limiting assumption holds,
r can be perfectly reconstructed from the discrete data{c[k]} byr(x) =∞
Xk=−∞esinc c)(x),c[k]sinc(x/∆ − k) =: (H

p

Q∆

r

c[k]

S∆

Q−1
∆

e sincHp̃

Fig. 2: Density estimation by Widrow’s method
e−ds
w

F (s)

p
Se∆

c

↑L

H∆
L

K(z)

+
−
p̃

e

Fig. 3: Error system E
where sinc is the sinc function defined bysinc(x) :=sin πx.
πxThen, the original PDF p(x) can be obtained by the−1−1convolution p = r ∗ q∆
, where q∆
is the inverse systemof q∆ , whose Fourier transform is given by (see (3))
Q−1
∆ (jω) =jω
22ejω ∆ − e−jω ∆.

This procedure is shown in Fig. 2.

4

H ∞ Optimal Estimation

4.1 Signal space F L2
In practice, the band-limiting assumption p ∈ BL2
can be too restrictive because many important probablity density functions such as uniform, Erlangian, orRayleigh ones do not satisfy the assumption (see [6]).
We therefore introduce another signal space F L2 defined by
F L2 := {p ∈ L2 : p = F w, w ∈ L2 },where F is a stable and strictly causal system with arational transfer function F (s). Although this spacedoes not include band-limited functions (since the sincfunction does not have a rational transfer function), thisspace includes many functions such as the distributions
mentioned above [6].
The function F (s) is treated as a smoothing parameter[1] in our density estimation. For selecting F (s), see
sub-section 4.5 and section 5.
Density estimation via sampled-data H ∞optimization
Assuming that the original PDF p(x) is in F L2 , weconsider the reconstruction problem by the error system shown in Fig. 3. In this figure, an L2 signal w isfiltered by F (s). This F (s) is a presumed frequency
domain model for the original PDF p. Then, this p is
4.2

sampled by the generalized sampler Se∆ defined by (2)and a discrete PDF {c[k]} is obtained. To reconstruct
p from the discrete data {c[k]}, we adopt the interpolation technique in multirate signal processing [11]. First,
the discrete data {c[k]} is upsampled by ↑ L, which isdefined byc[k], n = Lk, k = 0, 1, . . .
((↑ L)c)[n] =0,otherwise
Then the upsampled signal (↑ L)c is filtered by K(z)which we design. The output of this filter is converted
to a continuous function by the zero-order hold H∆/L
with period ∆/L, and we obtain an estimation p̃ ofthe original PDF p. In the error system Fig. 3, theestimation p̃ is compared with d-shifted version p(x−d)
of the original PDF p(x). The objective is to attenuate
the continuous reconstruction error e with respect tothe H ∞ norm of the error system
E := {e−ds − H∆/L K(z)(↑ L)Se∆ }F (s).

(4)

The optimal filter K(z) is then obtained by modern
sampled-data control theory, and has the following advantages:
1. the optimal filter is always stable;
2. the design takes the inter-sample behavior into
account;
3. the system is robust against the uncertainty of F ;4. the optimal FIR (finite impulse response) filter
is also obtainable via Linear Matrix Inequalities(LMI) assuming that the filter K(z) is an FIR
filter.
The last advantage is the key idea of our density estimation to consider constraints for PDF estimation. We
discuss this in the next section.
4.3 Constraints for PDF estimation
Since p(x) is a PDF, it satisfiesp(x) ≥ 0,
Z ∞p(x)dx = 1.

(5)
(6)

−∞

We therefore make the estimated PDF p̃(x) so as tomeet the constraints (5) and (6). We here assume thatthe filter K(z) is an FIR filter, that is,K(z) =NX

an z −n .

(7)

n=0

For mathematical reason, we define an = 0 for n >N or n < 0. Then, the constraints (5) and (6) canbe described by an LMI and an LME (linear matrix
equality).
Lemma 1 The estimation p̃(x) ≥ 0 for any x ∈ R if
and only ifan ≥ 0,n = 0, 1, . . . , N.

Since an = 0 if n > N or n < 0, we have

Proof
Let {d[n]} be the output of the filter K(z), that is,
d := K(↑ L)c. Then p̃(x) is the output of the zero-orderhold H∆/L for d[n], and it is clear that p̃(x) ≥ 0 for anyx ∈ R if and only if d[n] ≥ 0 for any n ∈ Z. Therefore,
we show d[n] ≥ 0 for any n ∈ Z if and only if an ≥ 0
for n = 0, 1, . . . , N .
First assume an ≥ 0, n = 0, 1, . . . , N . Then, we have∞Xd[n] = {(↑ Lc) ∗ a)[n] =c[k]an−Lk .

(8)

k=−∞

Since an = 0 for n > N and n < 0, and since c[k] ≥ 0
for all k ∈ Z, we have d[n] ≥ 0, for all n ∈ Z.
Conversely, assume that there exists a number l ∈ Zsuch that al < 0. Consider a discrete PDF(1,
k = 0,
c[k] =(9)0,
k 6= 0.
Then we haved[L + l] =∞X

c[k]aL+l−Lk = c[0]al = al < 0.

k=−∞2Lemma 2 Assume that
∆

∞
X

Z ∞

c[k] = 1.

(10)

if and only if
N
X

(11)

n=0

Proof
First, assume the equality (11). By using (8), we havep̃(x)dx =

−∞

=

∞
∆ X
d[n]
L n=−∞∞∞∆ X X
c[k]an−Lk
L n=−∞
k=−∞

=1−∆

∞
X

c[k]

k=−∞
∞
X

∞

X
∆
c[k]an−Lk
L n=−∞
k=−∞
( ∞
)
∞
X
∆ X
c[k]
=1+
an−Lk − L .
L
n=−∞
+

∞
X

k=−∞

an−Lk = L,

for all k ∈ Z. Hence (10) holds.
Conversely, assume the equality (12) does not hold.
Consider the discrete PDF {c[k]} is given by (9). Then
we have
( ∞
)
Z ∞
X
∆
p̃(x)dx = 1 +1=an − L .
L n=−∞−∞
Since ∆ > 0 and L > 0, this leads to contradiction. 2
By these lemmas, we can design the filter K(z) which
satisfies the constraints (5) and (6). Note that these
constraints are not satisfied in Widrow’s estimation.
4.4 Optimal FIR filter design
In this section, we give a design formula for our density estimation. Let a state space realization of F (s)
be


A B
F (s) =
(s).
C 0

Z ∆

eAθ dθ0Z ∞eA(k∆−τ −∆/2) Bw(τ )dτ.

0

DefineF∆ (s) :=an = L.

(12)

n=−∞

=C

−∞

Z ∞

if and only ifk∆p̃(x)dx = 1

an = L,

n=0

Then by (2), we have
Z (k+1)∆c[k] =p (x + ∆/2) dx

k=−∞

Then,

N
X



A
C∆

B
0



(s),

C∆ := CZ ∆eAθ dθ,0and assume that the reconstruction delay d satisfies
d = (m + 1/2) ∆. Then we have a generalized plant
representation for our H ∞ design shown in Fig. 4. In
this figure, e′ (t) = e(t − ∆/2), and the delay is an allpass function, the H ∞ norm will be the same.
For the sampled-data system shown in Fig. 4 can be
reduced to a discrete-time system via H ∞ discretization or fast-sample fast-hold (FSFH) method [2, 3, 5].
Assuming that the filter K(z) is an FIR filter givenin (7),we have the following discrete-time error system
Ed (z) for the sampled-data system (4) [3],Ed (z) = T1 (z) + W (a)T2 (z),where a := [a0 , a1 , . . . , aN −1 ] is the vector of the FIRcoefficients, and W (a) := a0 I a1 I . . . aN −1 I .
Thanks to the a-affine structure of Ed (z), a state space
representation of Ed (z) is given by


Ae
Be
Ed (z) =(z),
Ce (a) De (a)

e′

w


F (s)e−m∆sF∆ (s)


−1
0

12000

10000

8000

6000

4000

2000

H∆/L

S∆

↑L

0
−8

−6

−4

Fig. 4: Generalized plant for H ∞ optimal density estimation
where Ae and Be are constant matrices and Ce (a) andDe (a) affinely depend on the design parameter a. By
using the Kalman-Yakubovic-Popov lemma [4], our optimization can be described as follows. Find the optimal FIR coefficient a which minimizes γ > 0 subject
toP = P ⊤ ≻ 0,A⊤
e P Ae − P Be⊤ P Ae
Ce (a)A⊤e P Be−γI + Be⊤ P Be
De (a)

V (a)  0,

0

2

4

6

8

Fig. 5: Histogram of samples

K(z)

5



−2


Ce (a)⊤De (a)⊤  ≺ 0,−γI

trV (a) = L,

where V (a) = diag{a1 , a2 , . . . , aN −1 }. The inequality
V (a)  0 is for the constraint (5) by Lemma 1, and theequality trV (a) = L is for (6) by Lemma 2. These are
LMI’s and an LME, and the optimal coefficient a can beobtained effectively by using an optimization softwares
such as MATLAB. If the up-sampling ratio L is large(e.g., L = 32 or 64) and the number N of K(z) is also
large, the computation will be hard. In such a case, acutting plane method can be applied. See [12].
4.5 Selection of parameters ∆ and F (s)We have assumed that the quantization step size ∆and the smoothing parameter F (s) are given. However,in estimating a PDF, we have to select appropriate ∆and F (s) for good estimation. An optimal method forselecting ∆ is proposed in [9]. We can use this method
for selecting ∆. On the other hand, F (s) is typicallychosen as a low-pass filter with a cutoff frequency lessthan the Nyquist frequency π/∆. Although the cutoff
frequency can be arbitrarily chosen in theory, the reconstructed PDF will be very cheap if the cutoff frequency
is near or greater than π/∆. So the cutoff frequency
should be less than π/∆. A method for selecting F (s)is to estimate the cutoff frequency by interpolating andextrapolating the discrete Fourier transform of {c[k]}.

Numerical Example

We here show an example of density estimation. We
set the original PDF p(x) to be the uniform distribution in the interval [−5, 5]. Note that since the support of p(x) is [−5, 5] (compact), p(x) does not satisfy the band-limiting assumption of Widrow’s theorem.
The number of samples drawn from the PDF is 105 .
Then the samples are uniformly quantized with step
size ∆ = 1. The histogram of the quantized samples isshown in Fig. 5.
For our design, we set the design parameters as follows. The transfer function (smoothing parameter)
F (s) is set as
ω1 ω2, ω1 := π/10, ω2 := π/20.
F (s) =
(s + ω1 )(s + ω2 )The up-sampling ratio L is 4. The number of taps is
N = 16.
Fig. 6 shows the frequency response of H ∞ optimalfilter K(z), and Fig. 7 shows its impulse response. By
Fig. 7 we can see that our filter satisfies P
the constraintan ≥ 0 in Lemma 1. The constraint
an = L in
Lemma 2 means that K(1) = L (i.e., the frequency
response of K(z) is L at ω = 0), and Fig. 6 shows that
this constraint is satisfied.
Fig. 8 shows the estimated PDFs by Widrow’s
method and the proposed one. It can be seen that the
PDF of Widrow’s estimation has negative values. On
the other hand, our PDF is always positive. Moreover,
our method well reconstructs the original PDF, while
Widrow’s estimation shows large ripples. The ripples
are caused by the high frequency components in theoriginal characteristic function (F p)(ω). which is notconsidered in Widrow’s theorem. By this example, we
can see the efficiency of our method.

6

Conclusion

In this paper, we have proposed a new nonparametricdensity estimation from quantized samples via sampleddata H ∞ optimization. Our estimation is formulated
by signal reconstruction with an up-sampler and a digital filter, which can be effectively used in digital systems

Magnitude (dB)

20

0.12

0
0.1
−20
0.08

−40
−60
0

0.1

0.2

0.3
0.4
0.5
0.6
0.7
Normalized Frequency (×π rad/sample)

0.8

0.9

1

0.04

100
Phase (degrees)

0.06

0
0.02
−100
−200

0

−300
0

0.1

0.2

0.3
0.4
0.5
0.6
0.7
Normalized Frequency (×π rad/sample)

0.8

0.9

1

Fig. 6: Frequency response of H ∞ optimal filter K(z)
0.9

−0.02
−8

−6

−4

−2

0

2

4

6

8

Fig. 8: Density estimation: original PDF (dots),
Widrow’s estimation (dash), and H ∞ optimal one(solid)0.8[5] M. Nagahara and Y. Yamamoto, A new design forsample-rate converters, Proc. of 39th IEEE Conf.
on Decision and Control, pp. 4296–4301, 2000.

0.7
0.6
0.5

[6] A. Papoulis and S. U. Pillai, Probability, Random
Variables and Stochastic Processes, Fourth Ed.,
Mc Graw Hill, 2002.

0.4
0.3
0.2
0.1
0

2

4

6

8

10

12

14

16

[7] R. Schreier and G. C. Temes, Understanding DeltaSigma Data Converters, Wiley Interscience, 2005.
[8] C. E. Shannon, Communication in the presence ofnoise, Proc. IRE, Vol. 37, No. 1, pp. 10–21, 1949.

Fig. 7: FIR coefficients of H ∞ optimal filter K(z)such as DSP’s (Digital Signal Processors). The filter
design is formulated by an H ∞ optimization, and the
optimal filter can be computed by LMI’s and an LME.
We have studied density estimation of only one dimensional PDF. Future work includes higher dimensional
case.

References[1] C. M. Bishop, Pattern Recognition and Machine
Learning, Springer, 2006.
[2] T. Chen and B. A. Francis, Optimal Sampled-data
Control Systems, Springer, 1995.
[3] H. Ishii and Y. Yamamoto, Sampled-data H ∞ andH 2 /H ∞ design of multirate D/A converter, Trans.
of the ISCIE, Vol. 11, No. 10, pp. 585–592, 1998
(in Japanese).
[4] T. Iwasaki, LMI and Control, Shoko-do, 1997 (inJapanese).

[9] H. Shimazaki and S. Shinomoto, A method for selecting the bin size of a time histogram, NeuralComputation, Vol. 19, pp. 1503-1527, 2007.
[10] M. Unser, Sampling — 50 years after Shannon,
Proc. IEEE, Vol. 88, No. 4, pp. 569–587, 2000.
[11] P. P. Vaidyanathan, Multirate Systems and Filter
Banks, Prentice Hall, Englewood Cliffs, 1993.
[12] Y. Wakasa, D. Yasufuku, M. Nagahara, and Y. Yamamoto, Sampled-data design of interpolators using the cutting-plane method, Proc. of 41st IEEE
Conf. on Decision and Control, pp. 2564–2569,
2002.
[13] B. Widrow, A study of rough amplitude quantization by means of nyquist sampling theory, IRE
Trans. on Circuit Theory, Vol. 3, No. 4, pp. 266 –
276, 1956.
[14] B. Widrow, I. Kollar, and M.-C. Liu, Statisticaltheory of quantization, IEEE Trans. on Instrumentation and Measurement, Vol. 45, No. 2, pp.
353 – 361, 1996.