#
# 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.
#
#-----

source('utils.R')

# Ex 16.5:
#
n = 5
ds = c( 0.5, -1, 2 )
pnorm( -3 - ds * sqrt(n) ) - pnorm( 3 - ds * sqrt(n) ) + 1


# Ex 16.6:
#
# Construct the needed data for this problem (I could not find this data with the book):
#
n = 5
k = 22

sample_no = c()
for( ii in 1:k ){
    sample_no = c( sample_no, rep( ii, n ) )
}
data = c( 12.2, 12.1, 13.3, 13.0, 13.0,
          12.4, 13.3, 12.8, 12.6, 12.9,
          12.9, 12.7, 14.2, 12.5, 12.9,
          13.2, 13.0, 13.0, 12.6, 13.9,
          12.8, 12.3, 12.2, 13.3, 12.0,
          13.9, 13.4, 13.1, 12.4, 13.2,
          12.2, 14.4, 12.4, 12.4, 12.5,
          12.6, 12.8, 13.5, 13.9, 13.1,
          14.6, 13.4, 12.2, 13.7, 12.5,
          12.8, 12.3, 12.6, 13.2, 12.8,
          12.6, 13.1, 12.7, 13.2, 12.3,
          13.5, 12.3, 12.8, 13.1, 12.9,
          13.4, 13.3, 12.0, 12.9, 13.1,
          13.5, 12.4, 13.0, 13.6, 13.4,
          12.3, 12.8, 13.0, 12.8, 13.5,
          12.6, 13.4, 12.1, 13.2, 13.3,
          12.1, 12.7, 13.4, 13.0, 13.9,
          13.0, 12.8, 13.0, 13.3, 13.1,
          12.4, 13.2, 13.0, 14.0, 13.1,
          12.7, 12.4, 12.4, 13.9, 12.8,
          12.6, 12.8, 12.7, 13.4, 13.0,
          12.7, 13.4, 12.1, 13.2, 13.3 )

DF = data.frame( sample_no=sample_no, Moisture.Content=data )
res = control_limit_utils( DF$Moisture.Content, DF$sample_no )



# Make a control chart:
#
mu_0 = 13
sigma = 0.6

LCL = mu_0 - 3*sigma/sqrt(n)
UCL = mu_0 + 3*sigma/sqrt(n)

min_y = min( c( res$xbar_i_dot, LCL, UCL ) )
max_y = max( c( res$xbar_i_dot, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_6.eps", onefile=FALSE, horizontal=FALSE)
plot( res$xbar_i_dot, ylim=c(min_y, max_y), ylab='xbar_i' )
abline( h=mu_0, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid()
#dev.off()


# Ex 16.7:
#
LCL = res$x_bar_bar - 3 * res$sigma_est_from_std/sqrt(n)
UCL = res$x_bar_bar + 3 * res$sigma_est_from_std/sqrt(n)

min_y = min( c( res$xbar_i_dot, LCL, UCL ) )
max_y = max( c( res$xbar_i_dot, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_7.eps", onefile=FALSE, horizontal=FALSE)
plot( res$xbar_i_dot, ylim=c(min_y, max_y), ylab='xbar_i' )
abline( h=res$x_bar_bar, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid()
#dev.off()


# Ex 16.8:
#
LCL = res$x_bar_bar - 3 * res$sigma_est_from_range/sqrt(n)
UCL = res$x_bar_bar + 3 * res$sigma_est_from_range/sqrt(n)

min_y = min( c( res$xbar_i_dot, LCL, UCL ) )
max_y = max( c( res$xbar_i_dot, LCL, UCL ) )

#postsoript("../../WriteUp/Graphics/Chapter16/ex_8.eps", onefile=FALSE, horizontal=FALSE)
plot( res$xbar_i_dot, ylim=c(min_y, max_y), ylab='xbar_i' )
abline( h=res$x_bar_bar, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid()
#dev.off()


# Ex 16.9:
#
# Construct the needed data for this problem (I could not find this data with the book):
#
n = 6
xbar = c( 95.47, 97.38, 96.85, 96.64, 96.87, 96.52, 96.08, 96.48,
          96.63, 96.50, 97.22, 96.55, 97.02, 95.55, 96.29, 96.80,
          96.01, 95.39, 96.58, 96.43, 97.06, 98.34, 96.42, 95.99 )
ss = c( 1.30, 0.88, 1.43, 1.59, 1.52, 1.27, 1.16, 0.79,
        1.48, 0.80, 1.42, 1.65, 1.28, 1.14, 1.37, 1.40,
        1.58, .98, 1.21, .75, 1.34, 1.60, 1.22, 1.18 )
day = 1:length(xbar)

DF = data.frame( day=day, xbar=xbar, ss=ss )
x_bar_bar = mean( DF$xbar )
sigma_est_from_std = mean( DF$ss ) / get_a_sub_n(n)

LCL = x_bar_bar - 3 * sigma_est_from_std/sqrt(n)
UCL = x_bar_bar + 3 * sigma_est_from_std/sqrt(n)

min_y = min( c( DF$xbar, LCL, UCL ) )
max_y = max( c( DF$xbar, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_9.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$xbar, ylim=c(min_y, max_y), ylab='xbar_i' )
abline( h=x_bar_bar, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid()
#dev.off()


# Ex 16.10:
#
# Remove sample 22 and recompute
#
DF = DF[-22,]

x_bar_bar = mean( DF$xbar )
sigma_est_from_std = mean( DF$ss ) / get_a_sub_n(n)

LCL = x_bar_bar - 3 * sigma_est_from_std/sqrt(n)
UCL = x_bar_bar + 3 * sigma_est_from_std/sqrt(n)

min_y = min( c( DF$xbar, LCL, UCL ) )
max_y = max( c( DF$xbar, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_10.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$xbar, ylim=c(min_y, max_y), ylab='xbar_i' )
abline( h=x_bar_bar, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid()
#dev.off()


# Ex 16.12:
#
# Follow 16.6 (assume we have all of the variables from there loaded)
#
# Make a control chart (of the type we need for these additional tests):
#
mu_0 = 13
sigma = 0.6

LCL = mu_0 - 3*sigma/sqrt(n)
UCL = mu_0 + 3*sigma/sqrt(n)

min_y = min( c( res$xbar_i_dot, LCL, UCL ) )
max_y = max( c( res$xbar_i_dot, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_12.eps", onefile=FALSE, horizontal=FALSE)
plot( res$xbar_i_dot, ylim=c(min_y, max_y), ylab='xbar_i' )
abline( h=mu_0, col='black' )
abline( h=c( LCL, UCL ), col='red' )
abline( h=c( mu_0 - 2*sigma/sqrt(n), mu_0 + 2*sigma/sqrt(n) ), col='blue' )
abline( h=c( mu_0 - sigma/sqrt(n), mu_0 + sigma/sqrt(n) ), col='green' )
grid()
#dev.off()


# Ex 16.13:
#
# Follow 16.6 (assume we have all of the variables from there loaded)
#
# Make a control chart (of the type we need for these additional tests):
#
LCL = res$x_bar_bar - 3*res$sigma_est_from_IQR/sqrt(n)
UCL = res$x_bar_bar + 3*res$sigma_est_from_IQR/sqrt(n)
print( "Robust control limits:")
print( c(LCL, UCL) )