Skip to content

Instantly share code, notes, and snippets.

@amerberg amerberg/weibull-sum.R
Last active Aug 29, 2015

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
You can’t perform that action at this time.