sids<-c(2466,3317,3260,2807,1984,3232,2438,3941,3742,2892,3005,2495,2013,2722,
2807,3062,1616,3374,3062,2551,2863,3318,3033,4423,2722,3005,2977,2013,2098,
2353,3572,2495,2608,3118,3232,3175,2013,2750,3459,2353,2637,2863,3515,3515,2807,
3374,4394,1503)
plot(density(sids,n=200),type='l')
title('Sids data, N = 48')

# Compute bootstrap distribution of sample medians
m<-1000
med.boot<-numeric(m)
for (i in 1:m) med.boot[i]<-median(sample(sids,replace=T))
hist(med.boot,nclass=25,prob=T)
summary(med.boot)

#this was run three times storing results in x1, x2, and x3
plot(density(x1,n=200),type='l')
lines(density(x2,n=200),lty=2)
lines(density(x3,n=200),lty=2)
title('3 Bootstrap Distributions for the Median')

m<-1000
mean.boot<-numeric(m)
for (i in 1:m) mean.boot[i]<-mean(sample(sids,replace=T))
summary(mean.boot)

#this was run 3 times, storing the result in y1, y2, y3
plot(density(y1,n=200),type='l')
lines(density(y2,n=200),lty=2)
lines(density(y3,n=200),lty=2)
title('3 Bootstrap Distributions for the Mean')

#just for fun, plot medians against means
for (i in 1:m) {x<-sample(sids,replace=T)
med.boot[i]<-median(x)
mean.boot[i]<-mean(x)}
plot(mean.boot,med.boot,main='Medians vs Means, Bootstrap Values')
lines(lowess(mean.boot,med.boot),lwd=3)

#compute approximate 95% CI based on percentiles of the bootstrap
distributions
quantile(med.boot,c(.025,.975))
quantile(mean.boot,c(.025,.975))