#
# 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.
#
# Computations used for the problems in this section of the text.
#
#-----

save_plots = F

# Ex 6.1 (EPage 240):
#
DF = read.csv( "../../Data/CH06/ex6-01.txt", header=TRUE, quote="'" )
median( DF$Strength..MPa. )
sd( DF$Strength..MPa. )
sum( DF$Strength..MPa. > 10 )
sd( DF$Strength..MPa. ) / mean( DF$Strength..MPa. )

# Ex 6.3 (EPage 240):
#
DF = read.csv( "../../Data/CH06/ex6-03.txt", header=TRUE, quote="'" )
m = mean( DF$Thickness )

median( DF$Thickness )

s = sd( DF$Thickness )
z_crit = qnorm( 0.9 )
m + s * z_crit

pnorm( ( 1.5 - m ) / s )

# Part (e):
#
library(boot)
my_boot_median = function(d,bi){
    median( d[bi,] )
}
boot( DF, my_boot_median, 100 )

# Ex 6.4 (EPage 241):
#
Xs = read.csv( "../../Data/CH06/ex6-01.txt", header=TRUE, quote="'" )
colnames(Xs) = c("Strength")
Ys = read.csv( "../../Data/CH06/ex6-04.txt", header=TRUE, quote="'" )

mean(Xs$Strength) - mean(Ys$Strength)

V = var(Xs$Strength) / length(Xs$Strength) + var(Ys$Strength) / length(Ys$Strength)

sd(Xs$Strength) / sd(Ys$Strength)

var(Xs$Strength) + var(Ys$Strength)

# Ex 6.5 (EPage 241):
#
DF = read.csv( "../../Data/CH06/ex6-05.txt", header=TRUE, quote="'" )
x_bar = mean( DF$Audited )
y_bar = mean( DF$Book.value )
d_bar = mean( DF$Error )
N = 5000
T = 1761300

c( N * x_bar, T - N * d_bar, T*x_bar / y_bar )


# Ex 6.6 (EPage 241):
#
data = c( 127.96, 210.07, 203.24, 108.91, 178.21,
          285.37, 100.85, 89.59, 185.36, 126.94,
          200.19, 66.24, 247.11, 299.87, 109.64,
          125.86, 114.79, 109.11, 330.33, 85.54,
          117.64, 302.74, 280.55, 145.11, 95.36,
          204.91, 311.13, 150.58, 262.09, 477.08, 94.33 )
log_data = log(data)

mean( log_data )
var( log_data )

exp( mean( log_data ) + 0.5 * var( log_data ) )


# Ex 6.7 (EPage 241):
#
data = c( 103, 156, 118, 89, 125, 147, 122, 109, 138, 99 )
10000 * mean( data )

sum( data > 100 )

median( data )


# Ex 6.8 (EPage 241):
#
q_hat = 0.15
(1-q_hat)^2


# Ex 6.9 (EPage 241):
#
n = 150
DF = read.csv( "../../Data/CH06/ex6-09.txt", header=TRUE, quote="'" )

# Compute the mean of the full sample from the counts at each scratch level:
#
x_bar = sum( DF$Scratches * DF$Frequency ) / n


# Ex 6.11
#
n1 = 200
n2 = 200
x1 = 127
x2 = 176
x1 / n1 - x2 / n2

p1 = x1/n1
p2 = x2/n2

sqrt( p1*(1-p1) / n1 + p2*(1-p2) / n2 ) 


# Ex 6.14
#
n = 5
Xs = c( 237, 375, 202, 525, 418 )

max( Xs ) - min( Xs ) + 1


# Ex 6.15
#
data = c( 16.88, 10.23, 4.59, 6.66, 13.68, 14.23, 19.87, 9.4, 6.51, 10.95 )
sum( data^2 ) / ( 2 * length(data) )