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

source('utils.R')

# Ex 16.39:
#

# Calculate the limits of some control charts:
#
n = 6
k = 26
sum_xbar = 10980
sum_s = 402
sum_r = 1074

# Control charts for the process location:
# 
x_bar_bar = sum_xbar / k

sigma_est_from_std = ( sum_s / k ) * ( 1 / 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)
print( c(LCL, UCL) )

sigma_est_from_range = ( sum_r / k ) * ( 1 / get_b_sub_n(n) )
LCL = x_bar_bar - 3*sigma_est_from_range/sqrt(n)
UCL = x_bar_bar + 3*sigma_est_from_range/sqrt(n)
print( c(LCL, UCL) )

# Control charts for the process variation:
# 
sbar = sum_s / k
LCL = max( sbar - 3 * sbar * sqrt( 1 - get_a_sub_n(n)^2 ) / get_a_sub_n(n), 0 ) 
UCL =      sbar + 3 * sbar * sqrt( 1 - get_a_sub_n(n)^2 ) / get_a_sub_n(n) 
print( c(LCL, UCL) )

rbar = sum_r / k
LCL = max( rbar - 3 * get_c_sub_n(n) * rbar / get_b_sub_n(n), 0 ) 
UCL = rbar + 3 * get_c_sub_n(n) * rbar / get_b_sub_n(n)
print( c(LCL, UCL) )


# Ex 16.40:
#
data = c( 8, 1, 7, 5, 2, 0, 2, 3, 4, 3, 1, 2, 5, 7, 3, 4, 6, 5, 2, 4, 0, 10, 2, 6 )
xbar = mean( data )
LCL = max( xbar - 3 * sqrt( xbar ), 0 )
UCL = xbar + 3 * sqrt( xbar )
print( c( LCL, UCL ) )

min_y = min( c( data, LCL, UCL ) )
max_y = max( c( data, LCL, UCL ) )

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


# Ex 16.41:
#
n = 3
k = 22
sample_no = c()
for( ii in 1:k ){
    sample_no = c( sample_no, rep( ii, n ) )
} 

data = c( 51.3, 51.7, 49.5,
          51.0, 50.0, 49.3,
          50.8, 51.1, 49.0,
          50.6, 51.1, 49.0,
          49.6, 50.5, 50.9,
          51.3, 52.0, 50.3,
          49.7, 50.5, 50.3,
          51.8, 50.3, 50.0,
          48.6, 50.5, 50.7,
          49.6, 49.8, 50.5,
          49.9, 50.7, 49.8,
          49.6, 48.4, 50.0,
          49.8, 51.2, 49.7,
          50.4, 49.9, 50.7,
          49.4, 49.5, 49.0,
          50.7, 49.0, 50.0,
          50.8, 49.5, 50.9,
          48.5, 50.3, 49.3,
          49.6, 50.6, 49.4,
          50.9, 49.4, 49.7,
          54.1, 49.8, 48.5,
          50.2, 49.6, 51.5 )

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

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_41.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.42 (A np-chart): Data from Example 16.6
#
n = 100
DF = read.csv( "../../Data/CH16/exp16-6.txt", header=TRUE, quote="'" )
pbar = mean( DF$pi )

LCL = max( n * pbar - 3 * sqrt( n * pbar * ( 1 - pbar ) ), 0 )
UCL =      n * pbar + 3 * sqrt( n * pbar * ( 1 - pbar ) )

min_y = min( c( n * DF$pi, LCL, UCL ) )
max_y = max( c( n * DF$pi, LCL, UCL ) )

#postscript("../../WriteUp/Graphics/Chapter16/ex_42.eps", onefile=FALSE, horizontal=FALSE)
plot( n * DF$pi, ylim=c(min_y, max_y), ylab='n p_i' )
abline( h=n * pbar, col='black' )
abline( h=c( LCL, UCL ), col='red' )
grid ()
#dev.off()


# Ex 16.43:
#
DF = read.csv( "../../Data/CH16/ex16-45.txt", header=TRUE, quote="'" )

x_bar_bar = sum( DF$ni * DF$xi ) / sum( DF$ni ) 
sbar = sqrt( sum( ( DF$ni - 1 ) * DF$si^2 ) / sum( (DF$ni - 1) ) )

# A control chart for variation would have:
#
n = 3 # some of the samples 
LCL = max( sbar - 3 * sbar * sqrt( 1 - get_a_sub_n(n)^2 ) / get_a_sub_n(n), 0 ) 
UCL =      sbar + 3 * sbar * sqrt( 1 - get_a_sub_n(n)^2 ) / get_a_sub_n(n) 
print( c(n, LCL, UCL) )
n = 4 # other of the samples 
LCL = max( sbar - 3 * sbar * sqrt( 1 - get_a_sub_n(n)^2 ) / get_a_sub_n(n), 0 ) 
UCL =      sbar + 3 * sbar * sqrt( 1 - get_a_sub_n(n)^2 ) / get_a_sub_n(n) 
print( c(n, LCL, UCL) )

# A control chart for process location would have:
#
n = 3
LCL = x_bar_bar - 3*sbar/( get_a_sub_n(n) * sqrt(n) )
UCL = x_bar_bar + 3*sbar/( get_a_sub_n(n) * sqrt(n) )
print( c(LCL, UCL) )
n = 4
LCL = x_bar_bar - 3*sbar/( get_a_sub_n(n) * sqrt(n) )
UCL = x_bar_bar + 3*sbar/( get_a_sub_n(n) * sqrt(n) )
print( c(LCL, UCL) )


if( FALSE ){ 
    # Ex 16.44: Using data from Example 16.9 (data from Example 16.8):
    #
    # Compute the W_t = alpha \bar{X}_t + (1-\alpha) W_{t-1}
    #
    # mu_0 \pm 3 \sigma_t
    #
    DF = read.csv( "../../Data/CH16/exp16-8.txt", header=TRUE, quote="'" )

    alpha = 0.8
    mu_0 = 40
}