weibull-sum.R

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.

weibull.gif