Density 6. Estimation 1. Cross Validation 2. Histogram 3. Kernel Density Estimation 4. Local Polynomials 5. Higher Dimensions 6. Mixture Models 7. Converting Density Estimation Into Regression
6.1 Cross Validation Suppose we observe X1 , . . . , Xn from an unknown density f . Our goal is to estimate f nonparametrically. Finding the best estimator fn in some sense is equivalent to nding the optimal smoothing parameter h. How to measure the performance of fn ? Risk/Integrated Mean Square Error(IMSE) R(fn , f ) = E (L(fn , f )) where L(fn , f ) = (fn (x) f (x))2 dx. One can nd an optimal estimator that minimizes the risk function: fn = argmin R(fn , f ),
b fn
but the risk function is unknown!
1 2
Estimation of the Risk function Use leave-one-out cross validation to estimate the risk function. One can express the loss function as a function of the smoothing parameter h L(h) = = (fn (x) f (x)) dx (fn (x))2 dx 2
J(h) 2
Finding an optimal estimator fn E(J(h)) = E(J(h)) = R(h) + c The optimal smoothing parameter h = argmin J(h)
h
fn (x)f (x)dx +
f 2 (x)dx
Nonparametric estimator f can be expressed as a function of h and the best estimator fn can be obtained by Plug-in the optimal smoothing parameter f = f h
Cross-validation estimator of the risk function (up to constant)
2
J(h) =
fn (x)dx
2 n
n
f(i) (Xi )
i=1
where f(i) is the density estimator obtained after removing it h observation
3 4
6.2 Histogram Histogram Estimator Without loss of generality, we assume that the support of f is [0,1]. Divide the support into m equally sized bins 1 B1 = 0, m Let h =
1 m,
Theorem 1. Suppose that f is absolutely continuous and (f )2 du < . Then R(fn , f ) = h2 12 (f )2 du + 1 + o(h2 ) + o nh 1 n .
The optimal bandwidth is h =
1 2 , B2 = , m m
Bj
, . . . , Bm
n i=1
m1 = ,1 m
1 n1/3
6 (f (u))2 du
1/3
.
pj =
f (x)dx and Yj =
I(Xi Bj )
With the optimal binwidth, R(fn , f ) where C =
3 2/3 4
The histogram estimator is dened by
n
C n2/3
fn (x) =
j=1
pj I(x Bj ) h
(f (u))2 du
1/3
.
where pj =
Yj n . 5 6
Proof. For any x, u Bj , f (u) = f (x) + (u x)f (x) + for some x between x and u. Hence, pj =
Bj
Therefore, the bias of fn is (u x)2 f (x). 2 b(x) E (fn (x)) f (x) = = 1 E (pj ) f (x) h
f (u)du f (x) + (u x)f (x) +
Bj
=
(u x)2 f (x) du 2 x + O(h3 ).
pj f (x) h 1 1 f (x)h + hf (x) h j x + O(h3 ) f (x) = h 2 1 = f (x) h j x + O(h2 ). 2
= f (x)h + hf (x) h j
1 2
By the mean value theorem, for some xj Bj , b2 (x)dx
Bj
=
Bj
(f (x))2 h j
1 2 1 2
2
x x
dx + O(h4 )
2
= (f (xj ))2
Bj
h j
dx + O(h4 )
h3 + O(h4 ). = (f (xj ))2 12
7 8
Hence
1 m
Therefore, b2 (x)dx =
0 j=1 m Bj
b2 (x)dx + O(h3 )
0
1
m
m
v(x)dx
=
j=1 Bj m
v(x)dx =
j=1 Bj m j=1
pj (1 pj ) dx nh2 p2 j h2 f 2 (xj )
j=1
=
j=1
(f (xj ))2 h 12
2 1 0
h3 x + O(h3 ) 12
=
1 nh2
pj
j=1 Bj m
1 nh2
Bj
= For the variance of fn : v(x) Var (fn (x)) =
(f (x))2 dx + o(h2 ).
= =
1 1 nh nh2 1 1 nh nh
p2 = j
j=1 1
1 1 nh nh
1 pj (1 pj ) Var (pj ) = . 2 h nh2
f 2 (x)dx + o(1)
0
=
1 +o nh
1 n
By the mean value theorem, for some xj Bj , pj =
Bj
f (x)dx hf (xj ).
9
10
Theorem 2. The cross-validation estimator of risk for the histogram is m (n + 1) 2 p2 . J(h) = h(n 1) h(n 1) j=1 j Theorem 3. Let m = m(n) be the number of bins in the histogram fn . Assume that m(n) and m(n) = log n/n 0 as n . Dene
2 2
6.3 Kernel Density Estimation Given a kernel K and a positive number h, called the bandwidth, the kernel density estimator is: 1 fn (x) = n
n i=1
1 K h
x Xi h
.
ln (x) = where c =
max
z/(2m) 2
fn (x) c, 0
m n.
,
un (x) =
fn (x) c, 0
The choice of kernel K is not crucial but the choice of bandwidth h is important. We assume that K satises K(x)dx = 1,
2 xK(x)dx = 0 and K
Then
Pr ln (x) E (fn ) un (x)x 1 .
x2 K(x) > 0.
11
12
Theorem 4. Let Rx = E (f (x) f (x))2 be the risk at point x and let R = Rx dx denote the integrated risk. Assume that f is absolutely continuous and that int(f (x))2 dx < . Also, assume that K satises K(x)dx = 1, Then, Rx R = = f (x) K 2 (x)dx 1 4 4 + O(n1 ) + O(h6 ), K hn (f (x))2 + n 4 nhn K 2 (x)dx 1 4 4 K hn (f (x))2 dx + + O(n1 ) + O(h6 ). n 4 nhn
2 xK(x)dx = 0 and K
1 Proof. Write Kh (x, X) = h K
xX h
and fn (x) =
1 n
n i=1
Kh (x, Xi ).
E (fn (x)) = E (Kh (x, X)) xt 1 K = h h = =
f (t)dt
x2 K(x) > 0.
K(u)f (x hu)du K(u) f (x) huf (x) + h2 u2 f (x) + 2 du
1 = f (x) + h2 f (x) 2 Hence bias(fn (x)) = Similarly, Var (fn (x)) = f (x)
u2 K(u)du.
1 2 2 h f (x) + O(h4 ). 2 K n K 2 (x)dx + O(n1 ). nhn
14
13
Bandwidth selection The optimal smoothing bandwidth is h = c2 2 A(f )n c1
1/5
The Normal reference rule Assume that f is Normal, one can compute the h with the Gaussian kernel. h = 1.06n1/5 is estimated by = min{s, IQR/1.34} where s is the sample standard deviation and IQR is the interquartile range. The selected bandwidth is: hn = 1.06 . n1/5
where c1 = x2 K(x)dx, c2 = (K(x))2 dx and A(f ) = (f (x))2 dx. The only unknown quantity in h is A(f ) = (f (x))2 dx A(f ) = (f (x))2 dx = (f (4) (x))f (x)dx = E (f (4) ) where f (r) denote the r th derivative of f
15
16
Plug-in method Let r = E (f (r) ). To estimate 4 by using the Kernel method, one need to choose the optimal bandwidth which is a functional of 6 . 1. Estimate 8 with the bandwidth chosen the normal reference rule. 2. Estimate 6 with the bandwidth depending on 8 3. Estimate 4 with the bandwidth depending on 6 4. The selected bandwidth is h =
Cross-validation Cross-validation score function: J(h) = 2 f (x)dx n
2 n i=1
fi (Xi )
The selected bandwidth is h = argmin J(h)
h
c2 c2 4 n 1
1/5
Theorem 5 (Stones Theorem). Suppose f is bounded. Let fh denote the kernel estimator with bandwidth h and let h denote the bandwidth chosen by cross-validation. Then f (x) fh inf h
2
dx
2
a.s.
1.
f (x) fh
dx
17
18
Computation The cross-validation score function J(h) can be approximated by 1 J(h) = nh2
n n
Optimal Convergence rate From Theorem 4, the optimal convergence rate is O(n4/5 ) if the optimal bandwidth is used. Theorem 6. Let F be the set of all pdfs and let f (m) denote the mth derivative of f . Dene Fm (c) = For any estimator fn , sup
f Fm (c)
K
i=1 j=1
Xi X i h
+
2 K(0) + O nh
1 n2
where K (x) = K (2) (x) 2K(x) and K (2) K(z y)K(y)dy. For the computation of fn and J(h), use ; fast Fourier transform (FFT) binning strategy For the details, See Silverman (1981) and Wand and Jones (1995).
f F :
|f (m) (x)|2 dx c2 .
Ef
(fn (x) f (x))2 bn2m/(2m+1) .
where b > 0 is a universal constant that depends only on m and c.
19
20
Local log-likelihood 6.4 Local Polynomials Smoothed log-likelihood Recall that the nonparametric MLE of the pdf is p = (p1 , . . . , pn ) argmin L(p)
p
The kernel density estimator: fK (x) = argmax Lx The (p).
p
log-likelihood function of f :
n
L(f ) =
i=1
log f (Xi ) n
f (u)du 1 .
where L(p) =
n
n
log pi n
i=1 i=1
pi 1 .
The local log-likelihood at target value x is
n
A smoothed version of the log-likelihood at x (up to constant) is
n
Lx (f ) =
i=1
K
x Xi h
log f (Xi ) n
K
xu h
f (u)du.
Lx (p) =
i=1
K
x Xi h
n
log pi n
i=1
K
x Xi h
pi .
21
22
Local likelihood density estimator Approximate log f (u) with a polynomial
p
6.5 Higher Dimensions Suppose Xi = (Xi1 , . . . , Xid ). The multivariate kernel estimator is fn (x) = 1 n
n i=1
Px (a, u) =
j=0
aj (x u)j . j!
The local polynomial log-likelihood is
n
1 K(H1 (x X)). |H|
Lx (a) =
i=1
K
x Xi h
Px (a, Xi )n
K
u Xi h
ePx (a,u) du.
The local likelihood density estimator is
a a fn (x) = ePx (b,x) = eb0 ,
where a = (a0 , . . . , ap )T argmax Lx (a)
a
We often assume a simple form of the bandwidth matrix or kernel. For example, H =diag(h1 , . . . , hd ). Then d d xj Xij 1 K fn (x) = nh1 nd hj
i=1 j=1
23
24
The risk function; d 1 4 h4 R K j 4 j=1
2 fjj (x)dx + j=k
h2 h2 j k
The optimal bandwidth h = O(n1/(4+d) ). i
fjj fkk dx+
6.6 Mixture Models ( K (x)dx) . nh1 hd
2 d
The Normal Mixture Model
K
f (x, ) =
j=1
(xj )2 2 2 j
2 pj (x|j , j ),
In practice, one can choose the optimal bandwidth by cross-validation and often assume a simple form of the bandwidth matrix. H = h I. The curse of dimensionality.
where
2 j (x|j , j )
=
1 e 2
.
2 2 Dene = (p1 , . . . , pK , 1 , . . . , K , 1 , . . . , K )
Given K, how to estimate ? EM algorithm
25
26
EM algorithm For K = 2, Dene latent variables Zi : 1 : Zi = 0 : L(y, z, ) =
i:zi =0 n
[E-step] Given , compute E (Zi |Y, ). E(Zi |Y = yi , ) wi = if Yi is from group 1 if Yi is from group 2 where
n
p(yi ) (y) + (1 p)(y)
[M-step] = argmax L(),
The log-likelihood function is
2 log[p(yi |1 , 1 )]
+
i:zi =1
log[(1
2 p)(yi |2 , 2 )]
L() =
i=1
(wi log p(yi ) + (1 wi ) log(1 p)(yi )).
=
i=1
2 2 log[p(yi |1 , 1 )] + (1 zi ) log[(1 p)(yi |2 , 2 )]
Then
n
p =
2 1
=
i=1 n i=1
wi , n
1 =
n i=1 wi yi , n i=1 wi 2 2 =
2 =
n i=1 (1 wi )yi , n i=1 (1 wi
wi (yi 1 )2 , n i=1 wi
28
n i=1 (1 wi )(yi n i=1 (1 wi )
1 )2
.
27
How to estimate the number of components k? AIC (Akaike Information Criterion): nd most predictive model AIC = L q where L is the loglikelihood and p is the number of parameters. BIC (Bayesian Information Criterion): nd the true model with high probability q log n BIC = L 2
The number of clusters is not always equal to the number of component Gaussian is sensitive to outliers Replace the normal density with t-distribution density.
2 If j = 2 > 0 is xed and K n, then MLE of the mixture 1 model approaches the kernel estimate where pj = n and j = xj .
29
30
source("ch6.r") library(MASS) library(sm) library(locfit) library(mclust) ### Read data faithful<-read.table("faithful.dat",header=T) eruptions <-faithful$eruptions n<-length(eruptions) ### Select the optimal number of bins by CV hist.h <- cv.hist.fun(eruptions)$mbest ### Select optimal bandwidth of kernel estimators by normal ### reference rule/Cross-validation/plug-in sigma.hat <-min(sd(eruptions), IQR(eruptions)/1.34) h.normal <-1.06*sigma.hat/n^(0.2)
31
h.cv <- ucv(eruptions) h.plugin <-width.SJ(eruptions) f1<-density(eruptions,width=h.normal) # normal reference f2<-density(eruptions,width=h.cv) # cross-validation f3<-density(eruptions,width=h.plugin) # plugin f4<-density(eruptions,method="mclust") # finite mixture model f5 <-locfit(~eruptions, alpha=c(0.1,0.8),flim=c(1,6)) f6 <-locfit(~eruptions, alpha=c(0.1,0.6),flim=c(1,6), link="ident") postscript("density.ps") par(mfrow=c(2,2)) truehist(eruptions, xlim=c(1,6), ymax=0.8) # histogram truehist(eruptions, nbins=hist.h, xlim=c(1,6), ymax=0.8) plot(f1,ylim=c(0,0.8));plot(f2,ylim=c(0,0.8)) plot(f3,ylim=c(0,0.8));plot(f4,ylim=c(0,0.8)) plot(f5,ylim=c(0,0.8),main="local linear") plot(f6,ylim=c(0,...