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 <- lprde(x1,w1,lb=round(min(x))-0.5,binwidth=1,freq=TRUE)
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)
legend(max(x),ymax,xjust=1,yjust=1,cex=.8,
legend=c("N(34.5,1.5)", "SJ", "lprde(x,w)"),
col = c(1,1,2),
lty = c(1,2,1),
lwd = c(1,2,1))
Run the code above in your browser using DataLab