#
# file: Assignment4_Rsolution.R
#
# Author: B. M. Bolstad
# Date: Nov 16, 2004
#
# Aim: Sophisticated Solution code
#      
#
# Instructions: Copy and paste.
#
#
#

##
##
## part I
##
##

nsims <- 10000
samp.size <- 10
mean <- 100
sd <- 15


Sim.data <- matrix(rnorm(nsims*samp.size,mean,sd),nsims, samp.size)
Sim.sds <- apply(Sim.data,1,sd)
Sim.vars <- apply(Sim.data,1,var)

hist(Sim.sds,prob=TRUE)
lines(density(Sim.sds))

plot(ecdf(Sim.sds))

hist(Sim.vars,prob=TRUE)
lines(density(Sim.vars))

plot(edf(Sim.var))


sum(Sim.sds <=10)/nsims
sum(Sim.vars >= 600)/nsims


##
##
## Part II
##
##


# A Utility function:
# figure out what the MLE estimates would be if we fit the normal
# distribution to the specified data and then plot the density curve

plot.mle.norm <- function(x,add=FALSE){

  mle.mean <- mean(x)
  mle.var <- sum((x-mean(x))^2/length(x))

  x.pts <- seq(min(x),max(x),length=100)
  y.pts <- 1/sqrt(2*pi*mle.var)*exp(-1/(2*mle.var)*(x.pts- mle.mean)^2)

  if (add){
    lines(x.pts,y.pts,lwd=2,col="blue")
  } else {
    plot(x.pts,y.pts,type="l")
  }

}

plot.mle.norm(rnorm(1000,sd=1))


nsamps <- 1000


par(mfrow=c(2,2))

results <- matrix(0,8,2)
colnames(results) <- c("Mean","SD")
rownames(results) <-  c(1,2,4,8,16,32,64,128)

store.means <- matrix(0,nsamps,8)
colnames(store.means) <- c(1,2,4,8,16,32,64,128)


i <- 1
for (samp.size in c(1,2,4,8,16,32,64,128)){

  Sim2.data <- matrix(runif(nsamps*samp.size),nsamps,samp.size)

  means <- apply(Sim2.data,1,mean)
  hist(means,xlim=c(0,1),prob=TRUE,main=paste("Mean of",samp.size,"uniform r.v."))

  # good correspondance between red and blue means normal fitting well
  # at least at a rough visual level
  plot.mle.norm(means,add=TRUE)
  lines(density(means),lwd=2,col="red")


  store.means[,i] <- means
  
  results[i,1] <- mean(means)
  results[i,2] <- sd(means)
  i <- i+1
}


results

par(mfrow=c(1,1))
plot(c(0,0),xlim=c(0,1),ylim=c(0,15),type="n",xlab="x",ylab="density",main="Smoothed histograms of means")
palette <- rainbow(8)
for (i in 1:8){

  lines(density(store.means[,i]),col=palette[i],lwd=2)
}

legend(0.8,14,c(1,2,4,8,16,32,64,128),col=palette,lty=1,lwd=2)