Law of the Iterated Logarithm

LawIteratedLogarithm.R

limsup <- function(x){
 z = rep(0,length(x))
for(i in 1:length(x)){
 z[i] = max(x[i:length(x)])
 }
 return(z)
}


set.seed(2)
LIL <- function(  iter = 50000,
                  mu = 0,
                  sigma = 1){
    require(TeachingDemos)

        x <- rnorm(iter, mean = mu, sd = sigma)
        n <- 10:iter
        m <- sigma*sqrt(2*log(log(n))/n)
        plot( n, cumsum(x)[n]/n,
              log = "x",
              type = "n",
              ylab = expression(paste("Sample Mean ", bar(italic(X[n])))),
              xlab = "Cumulative Trials",
              xlim = c(10,iter),
              ylim = mu + c(-1,1))
        lines( n, m, lty = 2)
        lines( n, -1*m, lty = 2)
        xx <- c(n,rev(n))
        yy <- mu + c(m, -rev(m))
        polygon(xx, yy, col=grey(0.95), border = NA)
        abline( h = mu, lty = 2, lwd = 2, col="blue")
        lines(n, cumsum(x)[n]/n)
        mtext(expression(mu), adj = 1, side=4,cex=2, at = mu, col="blue")
        text( iter^0.75, mu+m[round(iter^0.75)], pos=3,labels = expression(mu %+-% sqrt( frac(2*sigma^2~plain(log)~plain(log)~n, n))))


plot.refresh <- function(...){
        iter <- slider(no = 1)
        mu <- slider(no = 2)
        sigma <- slider(no = 3)
        x <- rnorm(iter, mean = mu, sd = sigma)
        n <- 10:iter
        m <- sigma*sqrt(2*log(log(n))/n)
        plot( n, cumsum(x)[n]/n,
              log = "x",
              type = "n",
              main = "Law of the Iterated Logarithm",
              ylab = expression(paste("Sample Mean ", bar(italic(X))[n])),
              xlab = expression(paste("sample size ", italic(n), ", (log scale)", sep="")),
              xlim = c(10,iter),
              ylim = 0.5*c(-1,1))
        lines( n, mu + m, lty = 2)
        lines( n, mu - m, lty = 2)
        xx <- c(n,rev(n))
        yy <- mu+c(m, -rev(m))
        polygon(xx, yy, col=grey(0.95), border = NA)
        abline( h = mu, lty = 2, lwd = 2, col="blue")
        lines(n, cumsum(x)[n]/n)
        mtext(expression(mu), adj = 1, side=4,cex=2, at = mu, col="blue")
        text( iter^0.75, mu+m[round(iter^0.75)], pos=3,labels = expression( mu %+-% sqrt( frac(2*sigma^2~plain(log)~plain(log)~n, n))))

    }

    slider( plot.refresh,
            c("Sample size n",
              "Population mean",
              "Population standard deviation"),
            c(1000, -0.4, 0.5),
            c(100000, 0.4, 1.5),
            c(1, 0.1, 0.1),
            c(50000, 0, 1),
            title = "LIL Demo")
}