# hand written percentile bootstrap code.  
#  absolutely no claims that this is efficient

# heron data set is hard-coded here. That avoids some R stuff

nboot <- 1000    # number of bootstrap replicates

# create a vector to contain the bootstrap ests of ros mean
bootmean <- rep(0, nboot)

nobs <- length(heron$Pb)   # number of obs in data set

# loop to generate and process each bootstrap sample
for (iboot in 1:nboot) {
  i <- sample(1:nobs, repl=T)   # sample integer 1..N with replacement

  # make bootstrap set of observations
  xb <- heron$Pb[i]
  # with corresponding censoring flags
  cb <- heron$Cens[i]

  # calculate ROS estimate of mean from bootstrap sample
  #  this use of ros assumes logN distribution
  # provide arguments to ros() to change that

  bootmean[iboot] <- mean(ros(xb, cb))
  
# if wanted to avoid samples with too high a %censored, replace above by:

#  if (mean(cb) > 0.8) {bootmean[i] <- NA}
#   else bootmean[i] <- mean(ros(xb,cb))

#  without the comment symbols (#)

# end the loop 
   }

# report the 0.025 and 0.975 quantiles, after removing any NA's 
print(quantile(bootmean,c(0.025,0.975),na.rm=T))