R/myboot2.R

In Nat-Geo/math4753: A selection of lab functions

#' My bootstrap program
#'
#' Runs a bootstrap simulation to create a confidence interval
#'
#' Given a vector containing the population, creates iter random samples and runs a bootstrap to create a confidence interval of (1-alpha)%. Also produces a histogram of the resulting sampling distribution.
#'
#' @param iter the number of samples to take for a bootstrap simulation, default 10000
#' @param x a vector containing the population
#' @param fun a character string with the name of the function to be examined, default is "mean"
#' @param alpha the parameter determining the confidence level of the interval to 1-alpha percent, default is 0.05
#' @param cx a character expansion factor applied to the text added to the histogram
#' @param ... a series of parameters passed to the hist() call that generates the histogram
#'
#' @return ci, the confidence interval
#' @return x, the input population
#' @return fun, the function being examined
#' @export
#'
#' @examples
#'
#' set.seed(20)
#' sam=rnorm(20,mean=10,sd=4)
#' myboot2(x = sam)
myboot2<-function(iter=10000,x,fun="mean",alpha=0.05,cx=1.5,...){  #Notice where the ... is repeated in the code
  n=length(x)   #sample size

  y=sample(x,n*iter,replace=TRUE)
  rs.mat=matrix(y,nr=n,nc=iter,byrow=TRUE)
  xstat=apply(rs.mat,2,fun) # xstat is a vector and will have iter values in it
  ci=quantile(xstat,c(alpha/2,1-alpha/2))# Nice way to form a confidence interval
  # A histogram follows
  # The object para will contain the parameters used to make the histogram
  para=hist(xstat,freq=FALSE,las=1,
            main=paste("Histogram of Bootstrap sample statistics","\n","alpha=",alpha," iter=",iter,sep=""),
            ...)

  #mat will be a matrix that contains the data, this is done so that I can use apply()
  mat=matrix(x,nr=length(x),nc=1,byrow=TRUE)

  #pte is the point estimate
  #This uses whatever fun is
  pte=apply(mat,2,fun)
  abline(v=pte,lwd=3,col="Black")# Vertical line
  segments(ci[1],0,ci[2],0,lwd=4)      #Make the segment for the ci
  text(ci[1],0,paste("(",round(ci[1],2),sep=""),col="Red",cex=cx)
  text(ci[2],0,paste(round(ci[2],2),")",sep=""),col="Red",cex=cx)

  # plot the point estimate 1/2 way up the density
  text(pte,max(para$density)/2,round(pte,2),cex=cx)

  invisible(list(ci=ci,fun=fun,x=x))# Some output to use if necessary
}