#
# 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.19:
#
n = 100
k = 25
pbar = ( 578 / n ) / k
LCL = max( pbar - 3 * sqrt( pbar*(1-pbar) / n ), 0 )
UCL = pbar + 3 * sqrt( pbar*(1-pbar) / n )
print( c( LCL, UCL ) )

p_highest = 39/n
p_lowest = 13/n
print( c(p_highest, p_lowest) )


# Ex 16.20:
#
k = 30
n = 200
n_defective = c( 10, 18, 24, 17, 37, 19, 7, 25, 11, 24,
                 29, 15, 16, 21, 18, 17, 15, 22, 12, 20,
                 17, 18, 12, 24, 30, 16, 11, 20, 14, 28 )
pbar = ( sum( n_defective ) / n ) / k
LCL = max( pbar - 3 * sqrt( pbar*(1-pbar) / n ), 0 )
UCL = pbar + 3 * sqrt( pbar*(1-pbar) / n )
print( c( LCL, UCL ) )

ps = n_defective / n

# Plot the P chart:
#
min_y = min( c( ps, LCL, UCL ) )
max_y = max( c( ps, LCL, UCL ) )

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


# Ex 16.22: (have the variables from Ex 16.20 in memory)
#
ys = asin( sqrt( ps ) )
LCL = max( mean( ys ) - 3 * sqrt( 1 / (4*n) ), 0 )
UCL = mean( ys ) + 3 * sqrt( 1 / (4*n) )
print( c( LCL, UCL ) ) 


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


# Ex 16.25:
#
area_examined = c( .8, .6, .8, .8, 1.0, 1.0, .8, 1.0, .6, .6, .6, .8, .8, 1.0, 1.0, 1.0, .8, .6, .6, .6 )
num_blemishes = c( 3, 2, 3, 2, 5, 5, 10, 12, 4, 2, 1, 3, 5, 4, 6, 12, 3,3,5,1 )
us = num_blemishes / area_examined
ubar = mean( us )

# The control limits now depend on the point i:
#
LCL = ubar - 3 * sqrt( ubar / num_blemishes )
UCL = ubar + 3 * sqrt( ubar / num_blemishes )

# Find any out-of-control points:
#
bad_pts = which( ( us < LCL ) | ( us > UCL ) )
print( bad_pts )

print( area_examined[bad_pts] )
print( num_blemishes[bad_pts] )


# Ex 16.26:
#
ys = 2*sqrt( n_defects )
ybar = mean( ys )

LCL = max( ybar - 3, 0 )
UCL = ybar + 3
print( c( LCL, UCL ) )