Code for Drunk Article

Three R scripts were used to create the data and figures in the blog Extremely Deviant Yet Normal Drunks . These can be found on github at blog10118. Here I present the third script since it encompasses the two other scripts.


Creating Figures 3 and 4

# Rayleigh Flight of Drunk # Creating Fig 3 and 4 of blog101118 # Libraries library(ggplot2) library(data.table) # Parameters n <- 100 m <- 365 l <- 1000 sig <- rep(c(-1,1),2*n) s.dat <- data.table( k = numeric(), mk = numeric()) for (k in (1:l)){ y.dat <- data.table( j = numeric(), zj = numeric()) for (j in (1:m)){ u <- sample(sig,n) xn <- c(0, cumsum(u)) t <- seq(0, n, 1) d.dat <- data.table(cbind(t, xn)) zj <- d.dat[t==n,]$xn temp <- data.frame(cbind(j,zj)) y.dat <- rbind(y.dat, temp) } mk <- max(y.dat$zj) temp <- data.frame(cbind(k, mk)) s.dat <- rbind(s.dat, temp) } s.dat <- s.dat[order(-mk)] s.oep <- s.dat[,.N, by=mk] s.oep$ep <- cumsum(s.oep$N)/l write.csv(s.dat, "Rayleigh_AnnualMaxDisp_l1000.csv", row.names=FALSE) write.csv(s.oep, "Rayleigh_OEP_l1000.csv", row.names=FALSE) g <- ggplot() gp <- g + geom_point(aes(x=k, y=mk), data=s.dat) + xlab(expression(paste("Simulation year, ", italic(k)))) + ylab(expression(paste("Annual Maximum Displacement, ", italic(M[k])))) + theme_bw() + theme(text=element_text(family="Georgia", size=22)) ggsave(filename="AnnualDisp_Rayleigh_l1000.jpeg", plot=gp, width=8, height=6, dpi=200) g <- ggplot() gp <- g + geom_point(aes(x=mk, y=ep), data=s.oep, size=2) + geom_hline(yintercept=0.01, col="red") + geom_vline(xintercept=36.5, col="red") + xlab(expression(paste("Annual Maximum Displacement, ", italic(M[k])))) + ylab(expression(paste("Exceedance Probability, Pr {", italic(M[k])<=italic(y), "}"))) + scale_y_log10() + theme_bw() + theme(text=element_text(family="Georgia", size=22), panel.grid.minor=element_blank()) ggsave(filename="EP_Rayleigh_l1000.jpeg", plot=gp, width=8, height=6, dpi=200)