Skip to content

Instantly share code, notes, and snippets.

@jverzani
Created February 17, 2012 04:17
Show Gist options
  • Save jverzani/1850595 to your computer and use it in GitHub Desktop.
Save jverzani/1850595 to your computer and use it in GitHub Desktop.
Simluation code for manipulate
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")
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment