# Derived from https://github.com/tgouhier/climit
library(shiny)
library(openintro)
library(gridExtra)
library(BHH2)
seed = as.numeric(Sys.time())
shinyServer(function(input, output) {
rand_draw = function(n, p)
{
vals = NULL
vals = do.call(rbinom, list(n=n, size=1, prob=p))
return(vals)
}
rep_rand_draw = repeatable(rand_draw)
parent = reactive({
n = 1e5
return(rep_rand_draw(input$n, input$p))
})
samples = reactive({
pop = parent()
n = input$n
k = input$k
return(replicate(k, sample(pop, n, replace=TRUE)))
})
# plot 1
output$pop.dist = renderPlot({
popsize = 1000
counts = c(popsize*(1-input$p), popsize*input$p)/popsize
barplot(counts,
main= paste0("Population distribution: \n p = ", input$p),
#main = "Population distribution:",
names.arg=c(0,1),
ylab="Relative Frequency" ,ylim=c(0,1),
col=COL[1,2],
cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex.names = 1.5)
})
# plot 2
output$sample.dist = renderPlot({
par(mfrow=c(3,3))
x = samples()
par(mfrow=c(2,4))
for(i in 1:8){
counts <- table(x[,i])
barplot(counts, col = COL[7,2],
main = paste("Sample",i),
xlab = "", pch=19,
cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex.names = 1.5,
ylim = c(0,1.2*max(counts)))
box()
mean_samp = round(mean(x[,i]),2)
sd_samp = round(sd(x[,i]),2)
if(counts[1] > counts[2]){
legend("topright",
legend=bquote(atop(hat(p)[.(i)]==.(mean_samp))),
bty = "n", cex = 1.5, text.font = 3)
}
if(counts[1] <= counts[2]){
legend("topleft",
legend=bquote(atop(hat(p)[.(i)]==.(mean_samp))),
bty = "n", cex = 1.5, text.font = 3)
}
}
})
# text
output$num.samples = renderText({
k = input$k
paste0("... continuing to Sample ",k,".")
})
# plot 3
output$sampling.dist = renderPlot({
n = input$n
p = input$p
k = input$k
pop = parent()
ndist = colMeans(samples())
m_samp = round(mean(ndist),2)
sd_samp = round(sd(ndist),2)
ndens=density(ndist)
nhist=hist(ndist, plot=FALSE)
hist(ndist, main=paste("Sampling distribution:\n",
"Distribution of means of ", k,
" random samples, each\nconsisting of ", n,
" observations from the population", sep=""),
xlab="Sample means", freq=FALSE, ylim=c(0, max(ndens$y, nhist$density)),
col=COL[2,3], border = "white",
cex.main = 1.5, cex.axis = 1.5, cex.lab = 1.5)
legend_pos = ifelse(m_samp > 40, "topleft", "topright")
legend(legend_pos, inset = 0.025,
legend=bquote(atop("mean of " ~ hat(p)==.(m_samp),"sd of " ~ hat(p) ~ "(SE)" ==.(sd_samp))),
bty = "n", cex = 1.5, text.col = COL[2], text.font = 2)
lines(ndens, col=COL[2], lwd=3)
box()
})
# text
output$CLT.descr = renderText({
n = input$n ; p = input$p ; q = 1-p
pop = parent()
m_pop = p
n = input$n
se=round(sqrt(p*(1-p)/n),4)
paste("According to the Central Limit Theorem (CLT), the distribution of sample proportions
(the sampling distribution) should be nearly normal. The mean of
the sampling distribution is approximately equal to p (", m_pop, ")
and the standard error (the standard deviation of
sample proportions) should be approximately equal to the square root of probability of
success (p) times the probability of failure (1-p) divided by the sample size (sqrt(", p, "*", q,
"/",n, ") =", se,").")
})
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.