The bootstrap in the example is called a **non-parametric bootstrap**, or **case-resampling** (see [here][1] and [here][2], for example). 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. After drawing the new samples from your sample, you calculate the statistic of interest, this may be a mean, a median or regression coefficients. This is repeated $n$ times. In each iteration, some observations from your original sample are drawn multiple times while some observations may not be sampled at all. After the iterations, the mean, median and standard deviations of the $n$ bootstrap-estimates are calculated.

There are several methods for calculating confidence intervals based on the bootstrap samples. One very simple method is just calculating the 2.5 and 97.5% percentiles of the bootstrap samples (see code below).

Let's replicate the example from the website but using our own loop:

    #-----------------------------------------------------------------------------
    # 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]))
    
    original.estimates
    
    #-----------------------------------------------------------------------------
    # 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]))
      
    }
    
    #-----------------------------------------------------------------------------
    # The means, medians and SDs of the bootstrapped statistics
    #-----------------------------------------------------------------------------
    
    colMeans(store.matrix, na.rm=T)
    
    [1]  1.2997903  0.2526688 -1.5661686  0.2004957  0.9543834  0.2702361 -0.8996836  0.1788934
    [9]  2.0682618  4.0208757 -2.0969405  3.8660345
    
    apply(store.matrix,2,median, na.rm=T)
    
    [1]  1.3013079  0.2485721 -1.5571989  0.1986454  0.9252035  0.2688095 -0.9082379  0.1781161
    [9]  1.7718903  0.8269775 -1.6717102  0.6434859
    
    apply(store.matrix,2,sd, na.rm=T)
    
    [1]  0.39673804  0.03207790  0.26220013  0.01949437  0.48914523  0.02042504  0.22173921
    [8]  0.01667250  1.59102322 13.23433539  1.56310762 13.27525033
    
    #-----------------------------------------------------------------------------
    # The bootstrap bias is the difference between the mean bootstrap estimates
    # and the original estimates
    #-----------------------------------------------------------------------------
    
    colMeans(store.matrix, na.rm=T) - original.estimates
    
    [1] -0.0712912250 -0.0034460509 -0.0509238534  0.0049018955  0.0753405152  0.0009582592
    [7]  0.0856792712  0.0029512565  0.4654897865  3.1845709680 -0.4306843771  3.1870641941
    
    #-----------------------------------------------------------------------------
    # Quantile confidence intervals
    #-----------------------------------------------------------------------------
    
    apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T)
    
    [,1]      [,2]      [,3]      [,4]       [,5]      [,6]       [,7]      [,8]
    2.5%  0.5196003 0.1989776 -2.128998 0.1674420 0.03492612 0.2327157 -1.3066351 0.1449354
    97.5% 2.0604902 0.3229413 -1.091960 0.2417978 1.90245108 0.3137059 -0.4379761 0.2140118
    [,9]      [,10]     [,11]      [,12]
    2.5%  0.4414957  0.5811394 -8.434404  0.3363101
    97.5% 8.0471143 57.6416680 -1.115613 57.6061993

**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
    #-----------------------------------------------------------------------------
    
    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


   [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