#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
#-----

DF = read.csv("../../Data/chap_1_prob_11.csv")

# A histogram/density plot of the raw data (ignoring year):
#
#postscript("../../WriteUp/Graphics/Chapter1/ex_11_hist_plots.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,2))
hist( DF$MPG, main='Histogram of MPG' )
plot( density( DF$MPG ), main='KDE of MPG' )
par(mfrow=c(1,1))
#dev.off()

# A density plot one for each year:
#
mpg_range = range(DF$MPG)

#postscript("../../WriteUp/Graphics/Chapter1/ex_11_year_by_year_hists.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,3))
for( y in 1993:1995 ){
    mask = DF$Year == y
    plot( density( DF$MPG[mask] ), main=y, xlim=mpg_range )
    grid()
}
par(mfrow=c(1,1))
#dev.off()

# Bootstrap estimate of the precision of the median:
#
set.seed(1234)

# Perform bootstrap replicas:
#
n = length(DF$MPG)
B = 100

boot_medians = rep( NA, B )
for( bi in 1:B ){
    b = sample( DF$MPG, size=n, replace=TRUE )
    boot_medians[bi] = median(b)
}

m = median( DF$MPG )
print( m )
print( sd(boot_medians) )