Instantly share code, notes, and snippets.

Embed
What would you like to do?
library(matrixStats)
library(ggplot2)
library(animation)
#set seed for replicability
set.seed(1)
#generate 10^3 deviations, 10^6 replicates
w = replicate(10^6, rweibull(10^3, 0.9, 1))
#recall Weibull is subexponential when shape parameter is < 1.
#animate histogram of ratios of max to sum for first n
sums = rep(0, 10^6)
maxs = rep(0, 10^6)
increment = 10
saveGIF({
for (n in increment*(1:(1000/increment))) {
start = n-increment+1
end = n
#recalculate sum/max for each trial by adding `increment` new terms
sums = sums + colSums(w[start:end,])
maxs = colMaxs(rbind(w[start:end,], maxs))
print(qplot(maxs/sums, geom="histogram", main=sprintf("Histogram of Max/Sum (n=%d)",n),
xlab="ratio", xlim=c(0,1), ylim=c(0,3e5), binwidth=0.001))
}}, movie.name="weibull.gif", interval=0.2, ani.width=1000, ani.height=1000)
#Image shows ratio concentrated about 0 as number of terms increases to 1000.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment