setwd( "C:/Teaching/AMS-206-132/Winter-2017/NB10-JAGS" )

y <- c(  0,  2, 11, 18, 18, 24, 32, 33, 33, 36, 42, 44, 45, 45, 48, 48,
        48, 49, 50, 51, 55, 55, 56, 58, 59, 63, 65, 65, 68, 68, 69, 69,
        71, 72, 73, 74, 74, 74, 77, 77, 77, 78, 78, 79, 79, 79, 80, 80,
        80, 80, 81, 81, 81, 81, 82, 82, 82, 83, 83, 83, 84, 84, 85, 85,
        85, 86, 86, 86, 87, 87, 87, 88, 88, 89, 89, 90, 90, 90, 91, 91,
        92, 92, 93, 93, 95, 95, 95, 96, 96, 96, 97, 99, 99, 99, 99, 100,
       101, 101, 103, 104, 104 )

print( n <- length( y ) )

101

hist( y, breaks = 15, prob = T, xlab = 'y',
  ylab = 'Density', main = 'Final Exam Score' )

#####################################################################

qbinom( .025, 101, 0.5 )

41

qbinom( .975, 101, 0.5 ) 

60

# normal approximation to binomial:

round( n / 2 - qnorm( 0.975 ) * sqrt( n / 4 ) )

41

round( n / 2 + qnorm( 0.975 ) * sqrt( n / 4 ) )

60

c( y[ 41 ], y[ 60 ] )

# 77 83

median( y )

81

#####################################################################

plot.ecdf <- function(x, ..., ylab="Fn(x)", verticals = FALSE,
		      col.01line = "gray70", pch = 19)
{
    plot.stepfun(x, ..., ylab = ylab, verticals = verticals, pch = pch)
    abline(h = c(0,1), col = col.01line, lty = 2)
}

plot.ecdf( y, xlab = 'y', ylab = 'Empirical CDF of y',
  main = '' )

#####################################################################

M <- 10000

seed <- 1

set.seed( seed )

m.star <- rep( NA, M )

for ( i in 1:M ) {

  y.star <- sample( y, replace = T )

  m.star[ i ] <- median( y.star )

}

hist( m.star, prob = T, xlab = 'm', ylab = 'Density',
  main = 'Bootstrap Distribution of the Sample Median' )

quantile( m.star, probs = c( 0.025, 0.975 ) )

# 2.5% 97.5% 
#   77    83

par( mfrow = c( 2, 1 ) )

hist( y, breaks = 15, prob = T, xlab = 'y (Test Score)',
  ylab = 'Density', main = 'Final Exam Score', 
  xlim = c( 0, 104 ) )

hist( m.star, prob = T, xlab = 'm', ylab = 'Density',
  main = 'Bootstrap Distribution of the Sample Median',
  xlim = c( 0, 104 ) )

#####################################################################