Skip to content

Instantly share code, notes, and snippets.

@dmarcelinobr
Last active December 18, 2015 14:59
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save dmarcelinobr/5800912 to your computer and use it in GitHub Desktop.
Save dmarcelinobr/5800912 to your computer and use it in GitHub Desktop.
Annotated version of myboot
myboot<-function(data, stat, nreps, hist = TRUE) {
# Where 'data' is the original data set to be paste into 'myboot' program.
# 'stat' is the function that will generate the desired statistic such as "standard errors", "confidence intervals" etc.
# 'nreps' is the number of repetitions we want in the simulation.
estimates<-get(stat)(data)# Compute the number of estimates needed:
len<-length(estimates)
# Make a container object matrix for store the bootstrap results:
container<-matrix(NA, ncol = len , nrow = nreps)
# The length of "estimates" is the number of coefficients we will estimate standard errors for.
nobs<-nrow(data)#Compute the number of observations to resample
for(i in 1:nreps) {
# Now let's start our bootstrapping resampling loop
posdraws<-ceiling(runif(nobs)*nobs)
# Draw nobs uniform integers draws from 1 to nobs
resample<-data[posdraws,]
# This will sample a random sample from the original data
container[i,]<-get(stat)(resample)
# This apply the function 'stat' to the new data table and store it as a row in 'container' object
}
### 1.2 ### This will compute standard deviation for each column
sds<-apply(container,2,sd)
if(hist==T) {
mfrow=c(1,1)
frame() # alias for plot.new()
if(len<= 3) par(mfrow=c(len,1))
if((len> 3)&(len<= 6)) par(mfrow=c(3,2))
for(j in 1:len) hist(container[,j],
main=paste("Estimates for ", names(estimates)[j]), xlab="")
}
print(rbind(estimates,sds))
return(list(estimation=container, sds=sds))
}
###@ 2 @###
mod<-function(griliches76){
model<-lm(lw~s+iq, data=griliches76)
coef<-model[[1]]
# This function generates a vector with the boostrap estimates.
#mod1<-function(griliches76)lm(lw~s+iq, data=griliches76)[[1]]
}
mod1.res<-myboot(griliches76, "mod", 10000, hist=T)
# This will run the bootstrap function on the estimates 10.000 times and plot a histogram for the found distribution.
@briatte
Copy link

briatte commented Jun 24, 2013

Here's a remix of your code. The resampling is done with sample, and the function will show a progress bar. I have used the foreach package to make the computation parallelizable.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment