Simluation code for manipulate

simex.R

library(manipulate)
svalue <- identity

runSim <- function(nSamples, sizeSample, thePop, theStat, plotSamps, plotDist, plotPop) {
  m <- svalue(nSamples)
  n <- svalue(sizeSample)
  pop <- svalue(thePop)
  stat <- svalue(theStat)

  
    
    
    ## get samples
    samples <- matrix(NA,nrow=m, ncol=n)
    res <- numeric(m)
    for(i in 1:m) {
      if(grepl("^rnorm", pop))
        x <- rnorm(n)
      else if(grepl("^rt", pop))
        x <- rt(n, df=3)
      else if(grepl("^rexp", pop))
        x <- rexp(n)
      else if(grepl("^rbinom", pop))
        x <- rbinom(n, 1, 1/2) 

      samples[i,] <- x
      ## locally mask function
      prop <- mean(x)
      mean <- mean(x)
      median <- median(x)
      sd <- sd(x)
      IQR <- IQR(x)
      max <- max(x)
      min <- min(x)
      range <- max - min
      
      out <- try(eval(parse(text = stat)), silent=TRUE)
      if(inherits(out,"try-error")) {
        cat("Error with the statistic:",out,"\n")
        return(FALSE)
      }
      
      res[i] <- out
    }
    
    
    ## plot population
    if(svalue(plotPop)) {
      n <- 1000
      x <- try(eval(parse(text=pop)), silent=TRUE)
      if(inherits(x,"try-error")) {
        cat("Error with the population:",x,"\n")
        return(FALSE)
      }
      
      ## make densityplot
      dens <- density(x)
      mx <- max(dens$y)
      with(dens, plot(x,y, type="l", bty="n", lwd=2, ylim=c(-mx/10,mx)))
      with(dens, lines(x,y*0, col=gray(.8))) # add line
      if(svalue(plotSamps)) {
        ## add on some samples of size n with stat
        mx <- max(dens$y)
        for(i in 1:10) {
          points(samples[i,],rep(i*mx/11,svalue(sizeSample)), col=gray(.6))
          points(res[i], i* mx/11, cex=2, col="red", pch=15)
        }
        title(sprintf("Density plot of population with %s samples", 10))
      } else {
        title("Densityplot of population")
      }
    }
    
    ## do we plot sampling distr of statistic?
    if(svalue(plotDist)) {
      ## boxplot -- or densityplot
      if(svalue(plotPop)) {
        ## add as boxplot
        mx <- max(dens$y)
        boxplot(res, at=0, pars=list(boxwex=.1), add=TRUE,
                border="red", horizontal=TRUE)
      } else {
        dens <- density(res)
        with(dens, plot(x,y,
                        main="density plot of sampling distribution",
                        type = "l", col="red", bty="n"))
        rug(res[1:100], col="red")
      }
    }
}  

manipulate(
           runSim(nSamples, sizeSample, thePop, theStat, plotSamps,plotDist, plotPop),
           nSamples=picker(list(5, 10, 100, 1000), initial="100", label="Sample size"),
           sizeSample=picker(list(2,4,8,16,32,64,256), initial="32",  label="No. of samples"),
           thePop=picker(
             "rnorm(n, mean=0, sd=1)",
             "rt(n, df = 3)",
             "rexp(n, rate=1)",
             "rbinom(n,1,prob=1/2)",
             label="Population"),
           theStat=picker("mean","median",
             "sd","IQR",
             "max","min","range", label="Statistic"),
           plotPop = checkbox(FALSE, label="Plot population"),
           plotSamps=checkbox(TRUE, label="Plot some samples"),
           plotDist=checkbox(TRUE, label="Plot sampling distribution"),
           redo=button("Redo")
           )