H 2 /H ∞ Approach to the Histogram Method for Density Estimation
M. Nagahara†1 , K. I. Sato‡2 , and Y. Yamamoto†3†

Graduate School of Informatics, Kyoto University, ‡ Graduate School of Economics, Kyoto University
Kyoto, 606-8501, Japan1nagahara@ieee.org, 2 KenjiSato@e-mail.jp, 3 yy@i.kyoto-u.ac.jp

Abstract: In this paper, we study nonparametric density estimation by the histogram method. Histogram is interpreted as
quantization, which decreases the amount of information. Then interpolation (or estimation) of the missing informationis needed. To achieve this, we introduce sampled-data H 2 /H ∞ optimization. We design the reconstruction system which
optimizes the worst case error between the original PDF and the estimation. The optimization is formulated by linear
matrix inequalities and equalities. Numerical examples are illustrated to show effectiveness of our method.
Keywords: Density estimation, sampled-data control, histogram method.

1. INTRODUCTION
Probability density functions (PDF’s, henceforth) or
probability distributions play a fundamental role in bothanalysis 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 the other hand, because theanalysis and design are usually executed on a digital processor, acquired data must be discretized in both time andvalue and this severely reduces the information. In thispaper, we propose a method that weakens the latter difficulty. That is, we are concerned about estimating a PDF
from quantized data. Our analysis assumes no specificfunction form for the PDF, and hence our method is supposed to be classified into nonparametric estimation. At
the same time, 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 [14]. Thesepapers study the estimation problem from quantized samples by applying Shannon’s sampling theorem [10], [11].
We refer to the Widrow’s result as Widrow’s quantization
theorem. This theorem assumes that the original characteristic function is band-limited up to the frequency π/∆,where ∆ is the step size of the uniform quantizer. Thisfrequency corresponds to the Nyquist frequency in Shannon’s sampling theorem. The assumption which is calledband-limiting condition is quite restrictive because many
random variables such as uniform, exponential, Erlangian, or gamma ones do not satisfy the assumption [8].
The second helpful idea has come from the control theory. In order to avoid the stringent assumption of Shannon’s theorem, sampled-data control theory proved usefulin [6]. The paper shows that the sampled-data H 2 /H ∞optimization provides a way to obtain a filter without the
band-limiting assumption. By altering the method of [6]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 are unable toobtain the exact PDF. Such a restriction is inevitable 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 aninterpolator [12]. Like the method of [6], 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. In section 4, we
first point out problems in Widrow’s estimation. Then
to solve these problems, we propose a new estimation
method by sampled-data H ∞ optimization. We discuss amulti-dimensional case in section 5. Numerical examples
are given in section 6. Section 7 concludes our results.
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 , Fw denotes the Fouriertransform of w. We denote by χI the indicator function
on a set I ⊂ R; χI (x) = 1 if x ∈ I and χI (x) = 0 ifx∈
/ I. For a matrix M , M ⊤ and trM are respectively thetranspose and the trace of M . For a symmetric matrix P ,P ≻ 0, P  0 and P ≺ 0 denote positive, nonnegative,
and respectively negative definite matrices. For a linear
time-invariant (continuous-time or discrete-time) systemwith state space matrices {A, B, C, D},




A B
A B
(s) or(z)
C D
C D
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 theacquired samples are stored and processed in digital sys-

tems, and hence they must be converted to digital data,that is quantized. We assume that this quantization is a
uniform one with the step size ∆ and the no-overload input range ((−M − 1)∆, (M + 1)∆) (see [9]), which isdefined byp

Q : R → A := {−M ∆, (−M + 1)∆, . . . , M ∆},



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

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

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

w

r

Q∆

c[k]

S∆

p̃

Q−1
∆

esinc
H

Fig. 1 Density estimation by Widrow’s method
e−ds
F

+e
−

p
Se∆

+

↑L

K

H∆/L

n

Fig. 2 Block diagram for the optimal filter K.

(1)Then, our problem is formulated as follows:
Problem 1: Given samples x
e1 , xe2 , . . . , xeN which are
i.i.d. (independent, identically distributed) and quantizedby the uniform quantizer defined by (1), estimate the original PDF p(x) from which the samples are drawn.
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 that ofband-limited functions [14]. We introduce this solution
in the next section.

3. WIDROW’S QUANTIZATIONTHEOREMIn 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) andquantized by Q can be interpreted as samples of a discrete PDF c[k]. This discrete PDF is obtained by areasampling [14], denoted by Se∆ p, that is,
Z (k+1/2)∆
p(x)dx =: (Se∆ p)[k].
(2)
c[k] =(k−1/2)∆This can be rewritten byZ ∞χ[−∆/2,∆/2) (k∆ − x)p(x)dxc[k] =−∞= (q∆ ∗ p)(k∆) =: S∆ (q∆ ∗ p)[k],q∆ (x) := χ[−∆/2,∆/2) (x).
Thus, quantization can be interpreted as convolution byq∆ and ideal sampling S∆ . Note that the Fourier transform (frequency response) of q∆ is given by2
1  jω 2
Q∆ (jω) =e ∆ − e−jω ∆ .
(3)jω
Based on this idea, Widrow et al. proposed a reconstruction scheme based on Shannon’s sampling theorem. To
see this, we define the space BL2 of band-limited functions:
BL2 := p ∈ L2 : (Fp)(ω) = 0, |ω| > π/∆ .

Then, Widrow et al. have shown a solution to Problem
1 as follows [14].
Theorem 1 (Widrow et al.) Assume that p ∈ BL2 . Thenp(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) =: (Hwhere sinc(x) := sin(πx)/(πx). Then, the original PDF−1p(x) can be obtained by the convolution p = r ∗ q∆,−1
where q∆ is the inverse system of q∆ . This procedure isshown in Fig. 1.

4. H 2 /H ∞ OPTIMAL ESTIMATIONIn practice, the band-limiting assumption p ∈ BL2
can be too restrictive because many important probability density functions such as uniform, exponential, Erlangian, or gamma ones do not satisfy the assumption (see[8]). 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 a rational transfer function F (s). Although this space does
not include band-limited functions (since the sinc function does not have a rational transfer function), this space
includes many functions such as the distributions mentioned above [8]. The function F (s) is treated as asmoothing parameter [1] in our density estimation.
Assuming that the original PDF p(x) is in F L2 definedby we consider the reconstruction problem by the error
system shown in Fig. 2. In this figure, an L2 signal w
is filtered by F (s). This F (s) is a presumed frequency

domain model for the original PDF p. Then, this p is
sampled by the generalized sampler Se∆ defined by(Se∆ p)[k] :=Z (k+1/2)∆

p(x)dx(k−1/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. First, the discrete data {c[k]} is upsampled by ↑ L, which is definedby
c[l], k = Ll, l = 0, 1, . . .
((↑ L)c)[k] =0,otherwise
Then the upsampled signal (↑ L)c is filtered by K(z)which we design. The output of this filter is converted toa continuous function by the zero-order hold H∆/L withperiod ∆/L, and we obtain an estimation p̃ of the original PDF p. In the error system Fig. 2, the estimation p̃ iscompared with d-shifted version p(x − d) of the original
PDF p(x). The objective is to attenuate the continuousreconstruction error e with respect to the H ∞ norm of
the error system
J1 (K) := k{e−ds − H∆/L K(↑ L)Se∆ }F (s)kk∞ . (4)At the same time, to reduce the quantization noise n fromthe output, we limit the H 2 norm given by
J2 (K) := kH∆/L K(↑ L)k2 .

(5)

Our objective is to find K which minimizes J1 (K) subject to J2 (K) < λ for given λ > 0.
The sampled-data systems in (4) and (5) can be reduced to discrete-time ones via H ∞ discretization or fastsample fast-hold (FSFH) method [2], [3], [6]. Assumingthat the filter K(z) is an FIR filter
K(z) =N−1
X

al z −l ,l=0

we have the following discrete-time error systems E1 (z)and E2 (z) for the sampled-data systems in (4) and (5) [3],E1 (z) = T11 (z) + W (a)T12 (z),E2 (z) = W (a)T22 (z),where a := [a0 , a1 , . . . , aN −1
 ] is the vector of the FIRcoefficients, and W (a) := a0 I a1 I . . . aN −1 I .
The systems T11 , T12 , and T22 are all finite-dimensional
discrete-time LTI systems. Thanks to the a-affine structure of E1 (z) and E2 (z), their state-space representations
are given by


Ai
Bi
Ei (z) =(z), i = 1, 2,
Ci (a) Di (a)where Ai and Bi are constant matrices and Ci (a) andDi (a) affinely depend on the design parameter a. By

using the Kalman-Yakubovic-Popov lemma [4], our optimization can be described as follows: given λ > 0, findthe optimal FIR coefficient a which minimizes γ > 0
subject to
P1 = P1⊤ ≻ 0, P2 = P2⊤ ≻ 0
 ⊤A1 P1 A1 − P1
A⊤
1 P1 B1⊤
 B1 P1 A1−γI + B1⊤ P1 B1C1 (a)
D1 (a) ⊤A2 P2 A2 − P2 C2 (a)⊤≺ 0,
C2 (a)−I ⊤B2 P2 B2 − λ2 D2 (a)⊤≺ 0,
D2 (a)−I
V (a)  0,

C1 (a)⊤
D1 (a)⊤  ≺ 0,−γI

trV (a) = L,

where V (a) = diag{a1 , a2 , . . . , aN −1 }. The inequality
V (a)  0 is for the positivity of the estimated PDF p̃(x),that is, p̃(x) ≥ 0. The equality trV (a) = L is for theunity of the whole integral of p̃(x), that is,
Z ∞p̃(x)dx = 1.
−∞

For these conditions on V (a), see [7]. These are LMI’sand an LME, and the optimal coefficient a can be obtained 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 [13].

5. MULTI-DIMENSIONAL PDF
ESTIMATION
By using the one-dimensional filter K(z) given in theprevious section, we can construct a higher dimensional
PDF estimation system. In the case of M -dimensional
PDF’s, we adopt separable filter [5], that is,K(z1 , z2 , . . . , zM ) = K1 (z1 )K2 (z2 ) · · · KM (zM ),where each Ki (zi ) is a one-dimensional filter. By thisapproach, the computational complexity increases exponentially with the dimension M . In real system, this canbe applied when the dimension is at most 2 or 3.

6. NUMERICAL EXAMPLESWe 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. For our design, we set the design parameters as follows. The transfer function (smoothing parameter) F (s) is set as
F (s) =

ω1 ω2,
(s + ω1 )(s + ω2 )

0.12

0.1

p(x,y)
0.6

0.08

0.5
0.4

0.06

0.3
0.2

0.04

0.1
0.02

0

30
25

60
50

20

0

y

40

15

30

10
−0.02
−10

−8

−6

−4

−2

0

2

4

6

8

10

Fig. 3 Density estimation: original PDF (thin line),
Widrow’s estimation (dots), and proposed estimation(thick line)
26
24

20

5

x

10
0

Fig. 5 2-dimensional density estimation
Higher dimensional case have been also considered. Design examples have been illustrated to show effectivenessof the proposed method.

22

REFERENCES

20

y

18
16
14
12
10
8
10

15

20

25

30

35

40

45

50

x

Fig. 4 2-dimensional samples and histograms
where ω1 := π/10 and ω2 := π/20. The up-sampling
ratio L is 4. The number of taps is N = 16. Fig. 3
shows the estimated PDF’s 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, Widrow’s estimationshows large ripples, while our method well reconstructs
the original PDF. The ripples are caused by the high frequency components in Ψ(ω) which is not considered in
Widrow’s theorem. By this example, we can see the effectiveness of our method.
An example of 2-dimensional density estimation is
also shown. Fig. 4 shows given 2-dimensional data. From
the data, we estimate the original PDF. The obtained estimation is shown in Fig. 5.

7. CONCLUSION
In this paper, we have proposed a new nonparametricdensity estimation from quantized samples via sampleddata H ∞ optimization. Our estimation is formulated bysignal reconstruction with an up-sampler and a digital filter, which can be effectively used in digital systems suchas DSP’s (Digital Signal Processors). The optimizationof the filter is formulated by an H 2 /H ∞ one, and the
optimal filter can be computed by LMI’s and an LME.

[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).
[5] W.-S. Lu and A. Andoniou, Two-Dimensional Digital Filters, Marcel Dekker, 1992.
[6] M. Nagahara and Y. Yamamoto, A new design forsample-rate converters, Proc. of 39th IEEE CDC,
pp. 4296–4301, 2000.
[7] M. Nagahara, K. I. Sato, and Y. Yamamoto, H ∞ optimal nonparametric density estimation from quantized samples, Proc. of the 40th ISCIE SSS, 2008.
[8] A. Papoulis and S. U. Pillai, Probability, Random
Variables and Stochastic Processes, Fourth Ed., Mc
Graw Hill, 2002.
[9] R. Schreier and G. C. Temes, Understanding DeltaSigma Data Converters, Wiley Interscience, 2005.
[10] C. E. Shannon, Communication in the presence ofnoise, Proc. IRE, Vol. 37, No. 1, pp. 10–21, 1949.
[11] M. Unser, Sampling — 50 years after Shannon,
Proc. IEEE, Vol. 88, No. 4, pp. 569–587, 2000.
[12] P. P. Vaidyanathan, Multirate Systems and Filter
Banks, Prentice Hall, Englewood Cliffs, 1993.
[13] Y. Wakasa, D. Yasufuku, M. Nagahara, and Y. Yamamoto, Sampled-data design of interpolators using the cutting-plane method, Proc. of 41st IEEE
CDC, pp. 2564–2569, 2002.
[14] B. Widrow, I. Kollar, and M.-C. Liu, Statistical theory of quantization, IEEE Trans. on Instrumentationand Measurement, Vol. 45, No. 2, pp. 353 – 361,
1996.