There are several "flavours" or forms of the bootstrap (e.g. non-parametric, parametric, residual resampling and many more). The bootstrap in the example is called a **non-parametric bootstrap**, or **case resampling** (see [here][1], [here][2], [here][3] and [here][6] for applications in regression). The basic idea is that you treat your sample as population and repeatedly draw new samples from it **with replacement**. All original observations have equal probability of being drawn into the new sample. Then you calculate and store the statistic(s) of interest, this may be the mean, the median or regression coefficients **using the newly drawn sample**. This is repeated $n$ times. In each iteration, some observations from your original sample are drawn multiple times while some observations may not be drawn at all. After $n$ iterations, you have $n$ stored bootstrap estimates of the statistic(s) of interest (e.g. if $n=1000$ and the statistic of interest is the mean, you have 1000 bootstrapped estimates of the mean). Lastly, summary statistics such as the mean, median and the standard deviation of the $n$ bootstrap-estimates are calculated.

Bootstrapping is often used for:

 1. Calculation of confidence intervals (and estimation of the standard errors)
 2. Estimation of the bias of the point estimates

There are [several methods][4] for calculating confidence intervals based on the bootstrap samples ([this paper][5] provides explanation and guidance). One very simple method for calculating a 95%-confidence interval is just calculating the empirical 2.5th and 97.5th percentiles of the bootstrap samples (this interval is called the *bootstrap percentile interval;* see code below). The simple percentile interval method is rarely used in practice as there are better methods, such as the bias-corrected and accelerated bootstrap (BCa). BCa intervals adjust for both bias and skewness in the bootstrap distribution.

The *bias* is simply estimated as the difference between the mean of the $n$ stored bootstrap samples and the original estimate(s).

Let's replicate the example from the website but using our own loop incorporating the ideas I've outlined above (drawing repeatedly with replacement):

    #-----------------------------------------------------------------------------
    # Load packages
    #-----------------------------------------------------------------------------
    
    require(ggplot2)
    require(pscl)
    require(MASS)
    require(boot)
    
    #-----------------------------------------------------------------------------
    # Load data
    #-----------------------------------------------------------------------------
    
    zinb <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
    zinb <- within(zinb, {
      nofish <- factor(nofish)
      livebait <- factor(livebait)
      camper <- factor(camper)
    })
    
    #-----------------------------------------------------------------------------
    # Calculate zero-inflated regression
    #-----------------------------------------------------------------------------
    
    m1 <- zeroinfl(count ~ child + camper | persons, data = zinb,
                   dist = "negbin", EM = TRUE)
    
    #-----------------------------------------------------------------------------
    # Store the original regression coefficients
    #-----------------------------------------------------------------------------
    
    original.estimates <- as.vector(t(do.call(rbind, coef(summary(m1)))[, 1:2]))
        
    #-----------------------------------------------------------------------------
    # Set the number of replications
    #-----------------------------------------------------------------------------
    
    n.sim <- 2000
    
    #-----------------------------------------------------------------------------
    # Set up a matrix to store the results
    #-----------------------------------------------------------------------------
    
    store.matrix <- matrix(NA, nrow=n.sim, ncol=12)
    
    #-----------------------------------------------------------------------------
    # The loop
    #-----------------------------------------------------------------------------

    set.seed(123)
    
    for(i in 1:n.sim) {
      
      #-----------------------------------------------------------------------------
      # Draw the observations WITH replacement
      #-----------------------------------------------------------------------------
      
      data.new <- zinb[sample(1:dim(zinb)[1], dim(zinb)[1], replace=TRUE),]
      
      #-----------------------------------------------------------------------------
      # Calculate the model with this "new" data
      #-----------------------------------------------------------------------------
      
      m <- zeroinfl(count ~ child + camper | persons,
                    data = data.new, dist = "negbin",
                    start = list(count = c(1.3711, -1.5152, 0.879),
                                 zero = c(1.6028, -1.6663)))
      
      #-----------------------------------------------------------------------------
      # Store the results
      #-----------------------------------------------------------------------------
      
      store.matrix[i, ] <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
      
    }
    
    
    #-----------------------------------------------------------------------------
    # Save the means, medians and SDs of the bootstrapped statistics
    #-----------------------------------------------------------------------------
    
    boot.means <- colMeans(store.matrix, na.rm=T)
    
    boot.medians <- apply(store.matrix,2,median, na.rm=T)
    
    boot.sds <- apply(store.matrix,2,sd, na.rm=T)
    
    #-----------------------------------------------------------------------------
    # The bootstrap bias is the difference between the mean bootstrap estimates
    # and the original estimates
    #-----------------------------------------------------------------------------
    
    boot.bias <- colMeans(store.matrix, na.rm=T) - original.estimates
    
    #-----------------------------------------------------------------------------
    # Basic bootstrap CIs based on the empirical quantiles
    #-----------------------------------------------------------------------------
    
    conf.mat <- matrix(apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T),
    ncol=2, byrow=TRUE)
    colnames(conf.mat) <- c("95%-CI Lower", "95%-CI Upper")
    
And here is our summary table:

    #-----------------------------------------------------------------------------
    # Set up summary data frame
    #-----------------------------------------------------------------------------
    
    summary.frame <- data.frame(mean=boot.means, median=boot.medians,
    sd=boot.sds, bias=boot.bias, "CI_lower"=conf.mat[,1], "CI_upper"=conf.mat[,2])
    
    summary.frame
    
          mean  median       sd       bias CI_lower CI_upper
    1   1.2998  1.3013  0.39674 -0.0712912  0.51960   2.0605
    2   0.2527  0.2486  0.03208 -0.0034461  0.19898   0.3229
    3  -1.5662 -1.5572  0.26220 -0.0509239 -2.12900  -1.0920
    4   0.2005  0.1986  0.01949  0.0049019  0.16744   0.2418
    5   0.9544  0.9252  0.48915  0.0753405  0.03493   1.9025
    6   0.2702  0.2688  0.02043  0.0009583  0.23272   0.3137
    7  -0.8997 -0.9082  0.22174  0.0856793 -1.30664  -0.4380
    8   0.1789  0.1781  0.01667  0.0029513  0.14494   0.2140
    9   2.0683  1.7719  1.59102  0.4654898  0.44150   8.0471
    10  4.0209  0.8270 13.23434  3.1845710  0.58114  57.6417
    11 -2.0969 -1.6717  1.56311 -0.4306844 -8.43440  -1.1156
    12  3.8660  0.6435 13.27525  3.1870642  0.33631  57.6062

**Some explanations**

 - The difference between the mean of the bootstrap estimates and the original estimates is what is called "bias" in the output of `boot`
 - What the output of `boot` calls "std. error" is the standard deviation of the bootstrapped estimates

Compare it with the output from `boot`:

    #-----------------------------------------------------------------------------
    # Compare with boot output and confidence intervals
    #-----------------------------------------------------------------------------
    
    set.seed(10)
    res <- boot(zinb, f, R = 2000, parallel = "snow", ncpus = 4)
    
    res
    
    Bootstrap Statistics :
           original       bias    std. error
    t1*   1.3710504 -0.076735010  0.39842905
    t2*   0.2561136 -0.003127401  0.03172301
    t3*  -1.5152609 -0.064110745  0.26554358
    t4*   0.1955916  0.005819378  0.01933571
    t5*   0.8790522  0.083866901  0.49476780
    t6*   0.2692734  0.001475496  0.01957823
    t7*  -0.9853566  0.083186595  0.22384444
    t8*   0.1759504  0.002507872  0.01648298
    t9*   1.6031354  0.482973831  1.58603356
    t10*  0.8365225  3.240981223 13.86307093
    t11* -1.6665917 -0.453059768  1.55143344
    t12*  0.6793077  3.247826469 13.90167954

    perc.cis <- matrix(NA, nrow=dim(res$t)[2], ncol=2)
    for( i in 1:dim(res$t)[2] ) {
      perc.cis[i,] <- boot.ci(res, conf=0.95, type="perc", index=i)$percent[4:5] 
    }
    colnames(perc.cis) <- c("95%-CI Lower", "95%-CI Upper")
    
    perc.cis 

          95%-CI Lower 95%-CI Upper
     [1,]      0.52240       2.1035
     [2,]      0.19984       0.3220
     [3,]     -2.12820      -1.1012
     [4,]      0.16754       0.2430
     [5,]      0.04817       1.9084
     [6,]      0.23401       0.3124
     [7,]     -1.29964      -0.4314
     [8,]      0.14517       0.2149
     [9,]      0.29993       8.0463
    [10,]      0.57248      56.6710
    [11,]     -8.64798      -1.1088
    [12,]      0.33048      56.6702

    #-----------------------------------------------------------------------------
    # Our summary table
    #-----------------------------------------------------------------------------

    summary.frame
      
          mean  median       sd       bias CI_lower CI_upper
    1   1.2998  1.3013  0.39674 -0.0712912  0.51960   2.0605
    2   0.2527  0.2486  0.03208 -0.0034461  0.19898   0.3229
    3  -1.5662 -1.5572  0.26220 -0.0509239 -2.12900  -1.0920
    4   0.2005  0.1986  0.01949  0.0049019  0.16744   0.2418
    5   0.9544  0.9252  0.48915  0.0753405  0.03493   1.9025
    6   0.2702  0.2688  0.02043  0.0009583  0.23272   0.3137
    7  -0.8997 -0.9082  0.22174  0.0856793 -1.30664  -0.4380
    8   0.1789  0.1781  0.01667  0.0029513  0.14494   0.2140
    9   2.0683  1.7719  1.59102  0.4654898  0.44150   8.0471
    10  4.0209  0.8270 13.23434  3.1845710  0.58114  57.6417
    11 -2.0969 -1.6717  1.56311 -0.4306844 -8.43440  -1.1156
    12  3.8660  0.6435 13.27525  3.1870642  0.33631  57.6062

Compare the "bias" columns and the "std. error" with the "sd" column of our own summary table. Our 95%-confidence intervals are very similar to the confidence intervals calculated by `boot.ci` using the percentile method (not all though: look at the lower limit of parameter with index 9).

  [1]: http://www.stat.cmu.edu/~cshalizi/402/lectures/08-bootstrap/lecture-08.pdf
  [2]: http://bcs.whfreeman.com/ips5e/content/cat_080/pdf/moore14.pdf
  [3]: http://en.wikipedia.org/wiki/Bootstrapping_(statistics)#Case_resampling
  [4]: http://en.wikipedia.org/wiki/Bootstrapping_(statistics)#Deriving_confidence_intervals_from_the_bootstrap_distribution
  [5]: http://onlinelibrary.wiley.com/doi/10.1002/(SICI)1097-0258(20000515)19:9%3C1141::AID-SIM479%3E3.0.CO;2-F/abstract;jsessionid=2969C4F4A320EA6AEE19776DE70D4A5B.d02t04
  [6]: http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Bootstrapping.pdf