#Blood serum level
serum=c(1,1,2,4,4,6,8,8,9,9,13,13,14,15,16,16,17,17,17,18,
20,20,22,22,22,23,23,24,24,24,24,24,24,25,25,25,27,27,28,29,29,29,30,31,31,
32,32,33,33,35,44,45,64,94)
serum=serum/10
hist(serum)
mean(serum)
sd(serum)/sqrt(length(serum))
mean(serum,trim=.25)
#A simple nonparametric bootstrap algorithm for a default 25% trimmed mean
trim.boot=function(x,B=1000,trimp=.25){
nx=length(x)
tm=NULL
for(i in 1:B){
xs=sample(x,nx,replace=TRUE)
tmean=mean(xs,trim=trimp)
tm[i]=tmean
}
bt=mean(tm)
bs=sd(tm)
out=list(bt,bs)
names(out)=c("Bootstrap estimate of 25% Trimmed Mean","Bootstrap Standard Error")
out
}

# Computing bootstrap CI's by hand
npcitrim=function(x,trim=.25,B=1000,alpha=.05){
bootsamp=NULL
for(i in 1:B){
bootsamp[i]=mean(sample(x,length(x),replace=T),trim=trim)
}
that=mean(x,trim=trim)
z0=qnorm(sum(bootsamp<=that)/B)
bcl=pnorm(2*z0+qnorm(alpha/2))
bcu=pnorm(2*z0-qnorm(alpha/2))
bl=quantile(bootsamp,p=c(alpha/2,1-alpha/2,bcl,bcu))
names(bl)=c("Percentile Lower","Percentile Upper",
 "BC Lower","BC Upper")
return(bl)
}
#
# This script accompanies the trimmed mean portion of the 
# bootstrap handout
trim=function(x,d,trim=0){
return(mean(x[d],trim=trim))
}
b=boot(serum,trim,R=1000,stype="i",trim=.25)
boot.ci(b,type=c("perc","basic","bca"))