Talk:2005 French riots/stats

rm(list=ls) bootstrap <- function(data, B, theta, ...) { n <- NROW(data)   # number of observed samples cols <- NCOL(data) # number of columns in data frame boot <- list    # data structure to be returned boot$B <- B 	boot$data <- data boot$theta.hat <- theta(data, ...) # generate and store random sample of indices ind <- sample(1:n, n*B, replace=T) boot$indices <- matrix(ind, ncol=n, byrow=T) # make 3-D array of bootstrap sample data boot$samples <- aperm(array( t(data[ind,]), c(cols,n,B) ), c(2,1,3)) # compute theta.star of each bootstrap data sample boot$replicates <- apply(boot$samples, 3, theta, ...) # compute bootstrap estimate of standard error if(!is.vector(boot$replicates)) { boot$sderror <- sqrt(apply(boot$replicates,1,var)) } else { boot$sderror <- sqrt(var(boot$replicates)) } 	# compute bootstrap estimate of bias p.stars <- aperm(apply(boot$indices, 1, tabulate, n)/n, c(2,1)) p.bar <- apply(p.stars, 2, mean) theta.bar <- theta(data, wt=p.bar) boot$bias <- mean(boot$replicates) - theta.bar class(boot) <- "bootstrap" # object-orientation return(boot) } theta.hat <- function(boot.data, wt=NULL) { data = data.frame(day=1:6,cars=c(40,315,596,897,1295,1408)) data$cars <- data$cars + boot.data[,2] model <- lm(cars ~ day, data, weights=wt) prediction <- predict(model,newdata=data.frame(day=7)) return(prediction) } riot.data <- data.frame(day=1:6,cars=c(40,315,596,897,1295,1408)) model <- lm(cars ~ day, riot.data) B <- 250 # number of bootstrap replications to do boot.data <- data.frame(index=1:6, residuals=model$residuals) # boot residuals boot <- bootstrap(boot.data, B, theta=theta.hat) # replicate predictions library(MASS) truehist(boot$replicates,          main="Possible number of vehicles that will be burned on Nov 7", 	  xlab="bootstrap replicates of linear extrapolation",  	  ylab="observed frequency",  	  sub=c(paste("standard error =", round(boot$sderror, digits=2),  	       ", bias = ", round(boot$bias,digits=2))), col="light grey")
 * 1) Boostrap Algorithm
 * 1) theta-hat is prediction of number of cars to be burned overnight
 * 1) rioting data from French Ministry of Interior (number of vehicles burned)
 * 1) fit linear model
 * 1) plot histogram of replicates