#
# 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.14:
#
n = 8 
ss = c( .204, .315, .096, .184, .230, .212, .322, .287,
        .145, .211, .053, .145, .272, .351, .159, .214,
        .388, .187, .150, .229, .276, .118, .091, .056 )

sbar = mean( ss )
an = get_a_sub_n( n )
delta_s = 3 * sbar * ( sqrt( 1 - an^2 )/an )
LCL = max( 0, sbar - delta_s )
UCL = sbar + delta_s
s_control_1imits = c( LCL, UCL )
print( s_control_1imits )

# Plot the S chart:
#
min_y = min( c( ss, LCL, UCL ) )
max_y = max( c( 53, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_14.eps", onefile=FALSE, horizontal=FALSE)
plot( ss, ylim=c(min_y, max_y), ylab='s_i' )
abline( h=sbar, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid ()
#dev.off()

# Remove the out-of-control point
#
bad_index = which( ss == 0.388 )
ss = ss[-bad_index] # delete this point and rerun the code above


# Ex 16.15:
#
n = 4
k = 30
sum_of_ranges = 85.2

rbar = sum_of_ranges / k
bn = get_b_sub_n( n )
cn = get_c_sub_n( n )
delta_r = 3 * rbar * ( cn/bn )
LCL = max( 0, rbar - delta_r )
UCL = rbar + delta_r
r_control_1imits = c( LCL, UCL )
print( r_control_1imits )

n = 8
sum_of_ranges = 106.2 # run the above code again 


# Ex 16.16: Load the data from Ex 16.6 (in the R file ex16_02.R):
#

# Plot the S chart:
#
min_y = min( c( res$ss, res$s_control_limits[1], res$s_control_limits[2] ) )
max_y = max( c( res$ss, res$s_control_limits[1], res$s_control_limits[2] ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_16_s_chart.eps", onefile=FALSE, horizontal=FALSE)
plot( res$ss, ylim=c(min_y, max_y), ylab='s_i' )
abline( h=res$sbar, col='black' )
abline( h=c( res$s_control_limits[1], res$s_control_limits[2] ), col='red' )
grid()
#dev.off()

# Plot the R chart:
#
min_y = min( c( res$ss, res$r_control_limits[1], res$r_control_limits[2] ) )
max_y = max( c( res$ss, res$r_control_limits[1], res$r_control_limits[2] ) )
#postscript("../../WriteUp/Graphics/Chapter16/ex_16_R_chart.eps", onefile=FALSE, horizontal=FALSE)
plot( res$srange, ylim=c(min_y, max_y), ylab='r_i' )
abline( h=res$rbar, col='black' )
abline( h=c( res$r_control_limits[1], res$r_control_limits[2] ), col='red' )
grid()
#dev.off()


# Ex 16.17: Load the data from Ex 16.9 (in the R file ex16_02.R):
#
sbar = mean( DF$ss )
an = get_a_sub_n( n ) # n=6 for this problem
delta_s = 3 * sbar * ( sqrt( 1 - an^2 )/an )
LCL = max( 0, sbar - delta_s )
UCL = sbar + delta_s
s_control_limits = c( LCL, UCL )
print( s_control_limits )

# Plot the S chart:
#
min_y = min( c( DF$ss, s_control_limits[1], s_control_limits[2] ) )
max_y = max( c( DF$ss, s_control_limits[1], s_control_limits[2] ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_17_S_chart.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$ss, ylim=c(min_y, max_y), ylab='s_i' )
abline( h=sbar, col='black' )
abline( h=c( s_control_limits[1], s_control_limits[2] ), col='red' )
grid()
#dev.off()


# Ex 16.18: Load the data from Ex 16.9 (in the R file ex16_02.R):
#
star = mean( DF$ss^2 )
chi_crit_L = qchisq( 0.001, n-1 )
LCL = star * chi_crit_L / (n-1)
chi_crit_U = qchisq( 1-0.001, n-1 )
UCL = star * chi_crit_U / (n-1)
print( c( LCL, UCL ) )