[R] density() with confidence intervals

[David Croll]

Thank you very much for your help, Greg!

Here's my ludicrous vision/attempt/whatever :o)


### Idea for making a density() with confidence interval ###

xx <- faithful$eruptions

xx.hist <- hist(xx, breaks="FD", freq=F)

# plot(xx.hist$mids, xx.hist$density)
# gives a rough "plot(density(xx))" output

# xx.hist$counts is later used to "repeat" observations
# for loess(), making the CI band narrower where there
# are more observations (or data in each histogram "bin")

x <- c() # preparing variables...
y <- c()

for (i in 1:length(xx.hist$mids)) {
	# going through each xx.hist$mids/xx.hist$density pair
	for (j in 1:xx.hist$counts[i]) {
	# ...and repeating observations according to xx.hist$counts
		x <- append(x, xx.hist$mids[i])
		y <- append(y, xx.hist$density[i])		
		}
	}

xx.loess <- loess(y ~ x, span=1)
xx.predict <- predict(xx.loess, se=T)

plot(x,y)
lines(x, xx.predict$fit, col="red")
lines(x, xx.predict$fit + 1.96 * xx.predict$s, col="red", lty=5)
lines(x, xx.predict$fit - 1.96 * xx.predict$s, col="red", lty=5)


# Kind regards, David

> Here is a simple approach that uses bootstrapping (this could
> probably be improved by using better bootstrap estimates and not
> ignoring the dependence between points):
> 
> xx <- faithful$eruptions
> 
> fit1 <- density(xx)
> 
> fit2 <- replicate(10000, { x <- sample(xx, replace=TRUE); density(x,
> from=min(fit1$x), to=max(fit1$x))$y } )
> 
> fit3 <- apply(fit2, 1, quantile, c(0.025,0.975) )
> 
> plot(fit1, ylim=range(fit3)) polygon( c(fit1$x, rev(fit1$x)),
> c(fit3[1,], rev(fit3[2,])), col='grey', border=F) lines(fit1)
> 
>