mu = 34.5; s=2.5; n = 10000
x = round(rnorm(n, mu, s),1)
x0 = seq(min(x)-s,max(x)+s, length=100)
f0 = dnorm(x0,mu, s); ymax <- max(f0*1.2)
xt = table(x); n = length(x)
x1 = as.numeric(names(xt))
w1 = as.numeric(xt)
est1 <- wkde(x1,w1, bandwidth='nrd0',freq=TRUE)
est2 <- wkde(x1,w1, bandwidth='nrd',freq=TRUE)
est3 <- wkde(x1,w1, bandwidth='amise',freq=TRUE)
est4 <- wkde(x1,w1, bandwidth='mise',freq=TRUE)
est5 <- wkde(x1,w1, bandwidth='lscv',freq=TRUE)
est6 <- wkde(x, bandwidth='lscv')
est0 = density(x1,bw="SJ",weights=w1/sum(w1));
plot(f0~x0, xlim=c(min(x),max(x)), ylim=c(0,ymax),
xlab="x", ylab="Density", type="l")
lines(est0, col=1, lty=2, lwd=2)
lines(est1, col=2)
lines(est2, col=3)
lines(est3, col=4)
lines(est4, col=5)
lines(est5, col=6)
lines(est5, col=6,lty=2)
legend(max(x),ymax,xjust=1,yjust=1,cex=.8,
legend=c("N(34.5,1.5)", "SJ", "nrd0",
"nrd","amise","mise","lscv(x,w)","lscv(raw)"),
col = c(1,1,2,3,4,5,6,6),
lty = c(1,2,1,1,1,1,1,2),
lwd=c(1,2,1,1,1,1,1,1))
Run the code above in your browser using DataLab