# boottest.r                 December 2007
 #
 # compare speeds for bootstrap
 #
 #
 # set up
 #
 #  Exercise 8.13 -- least 3/2's estimator
 #
 #  define function
  g3o2 <- function(x) {
 #                  define function to be minimized
    f3o2 <- function (mu)  sum((abs(x-mu)**(3/2))) 
 #                  derivative of function to be minimizes
    fp3o2 <- function (mu,u) { # derivative  
            fp3o2 <- sum( sign(u-mu)*sqrt(abs(u-mu)) ) }
 #  
 #  find root of derivative function
    that2 <- uniroot(f=fp3o2, lower=min(x), upper=max(x), u=x )
    g3o2 <- that2$root        }    # root is statistic
 #
 #  data
  v <- c(225,215,209,175,163,135,153,125)
 #
  print(g3o2(v))        # compute statistic on data
 # First method -- create matrix and use apply
 # data
  set.seed(5151917)
  Nboot <- 1000
  g3o2(v)               # compute statistic on data
 # 
  proc.time()       # read clock before execution
 #
 #       sample from v -- how many? length(v)*Nboot, WITH replacement
 #     matrix  so that each column is a bootstrap sample
  V <- matrix(sample(v,length(v)*Nboot,replace=T),length(v),Nboot)
 # compute statistic
  vps <- apply(V,2,g3o2)
  proc.time()       # read clock after execution
  var(vps)
  mean(vps)
 #
 # Second method -- using loop
  lv <- length(v)
  set.seed(5151917)                # same random numbers
  proc.time()       # read clock before execution
  vps <- rep(0,Nboot)              # initialize to zero
  for(i in(1:Nboot)){              # start loop
  vbs <- sample(v,lv,replace=T)
  vps[i]<-g3o2(vbs)           }    # end loop
  proc.time()       # read clock after execution
  mean(vps)
  var(vps)
 # done
  rm(list=ls())
 q()