[R] Bootstrap

[Matt Shotwell]

In order to apply the bootstrap, you must resample, uniformly at random
from the independent units of measurement in your data. Assuming that
these represent the rows of 'data', consider the following:

est <- function(y, x, obeta = c(1,1), verbose=FALSE) {
    n <- length(x)
    X <- cbind(rep(1, n), x)
    nbeta <- c(0,0)
    iter <- 0
    while(crossprod(obeta-nbeta)>10^(-12)) {
        nbeta <- obeta
        eta   <- X%*%nbeta
        mu    <- eta
        mu1   <- 1/eta
        W     <- diag(as.vector(mu1))
        Z     <- X%*%nbeta+(y-mu)
        XWX   <- t(X)%*%W%*%X
        XWZ   <- t(X)%*%W%*%Z
        Cov   <- solve(XWX)
        obeta <- Cov%*%XWZ
        iter  <- iter+1
        if(verbose)
            cat("Iteration #  and beta1= ",iter, nbeta, "\n")
    }
    return(nbeta[1,1])
}

boot <- function(data, reps) {
    n <- nrow(data)
    Nt <- vector('numeric', length=reps)
    for(Ncount in 1:reps) {
        #resample the rows of data
        bdata <- data[sample(1:n,n,replace=TRUE),]
        #recompute and store estimate
        Nt[Ncount] <- est(bdata[,1], bdata[,2])
    }
    return(Nt) 
}

stem(boot(data,1000),width=60)

  The decimal point is at the |

  -3 | 4
  -2 | 
  -1 | 2
  -0 | 888666666555444333333333332222222222211111111111
   0 | 000000000000001111111111111111111111111111222222+400
   1 | 000000000000000000000000000000011111111111111111+203
   2 | 000000000000000011111111222222222223333444444444+23
   3 | 000000001111111111222222222223344444445555556666
   4 | 11222233445556666789
   5 | 02334446677899
   6 | 1112334455778
   7 | 000011235568
   8 | 001799
   9 | 0259
  10 | 1446
  11 | 19
  12 | 48
  13 | 8
  14 | 024
  15 | 
  16 | 
  17 | 0788
  18 | 
  19 | 1

On Wed, 2011-07-20 at 18:09 -0400, Val wrote:
> Hi all,
> 
> I am facing difficulty on  how to use bootstrap sampling and
> below is my example of function.
> 
> Read a data , use some functions and  use iteration to find the solution(
> ie, convergence is reached).  I want to use bootstrap approach to do it
> several times (200 or 300 times) this whole process  and see the
> distribution of parameter of interest.
> 
> Below is a small example that resembles my problem. However,  I  found out
> all samples are the same. So I would appreciate your help on this case.
> 
> #**************************************
> rm(list=ls())
>  xx <- read.table(textConnection(" y x
>     11 5.16
>     11 4.04
>     14 3.85
>     19 5.68
>     4 1.26
>     23  7.89
>     15 4.25
>     17 3.94
>     7 2.35
>     17 4.74
>     14 5.49
>     11 4.12
>     17 5.92"), header=TRUE)
>     data <- as.matrix(xx)
>     closeAllconnections()
> 
> Nt <- NULL
> for (Ncount in 1:100)
>  {
>     y <- data[,1]
>     x <- data[,2]
>     n <- length(x)
> 
>     X <- cbind(rep(1,n),x)                 #covariate/design matrix
>     obeta<- c(1,1)                         #previous/starting values of beta
> 
>     nbeta <- c(0,0)                        #new beta
>     iter=0
> 
>   while(crossprod(obeta-nbeta)>10^(-12))
>    {
>     nbeta <- obeta
>     eta   <- X%*%nbeta
>     mu    <- eta
>     mu1   <- 1/eta
>     W     <- diag(as.vector(mu1))
>     Z     <- X%*%nbeta+(y-mu)
>     XWX   <- t(X)%*%W%*%X
>     XWZ   <- t(X)%*%W%*%Z
>     Cov   <- solve(XWX)
>     obeta <- Cov%*%XWZ
>     iter  <- iter+1
> 
>     cat("Iteration #  and beta1= ",iter, nbeta, "\n")
>     }
> 
>   Nt[Ncount] <- nbeta[1,1]
> }
> Nt
> summary(Nt)
> #**************e*****************************************
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.