R/main.R
In bolus123/PH1AND2XBAR: Phase I and Phase II Shewhart X-bar Chart
Defines functions
W.f
W.moment
pars.root.finding
d2.f
c4.f
ub.cons.f
PH1.corr.f
PH1.get.cc.mvt
PH1.joint.pdf.mvn.chisq
PH1.root.mvn.F
PH1.get.cc.mvn
PH1.rmvn.MC
PH1.joint.pdf.mvn.chisq.MC
PH1.root.mvn.F.MC
PH1.get.cc.mvn.MC
PH1.get.cc
PH1XBAR
PH1XBAR.data
PH2.inner.normal
PH2.CFAR.intgrand
PH2.CARL.intgrand
PH2.get.cc.uc
PH2.root.finding.EPC
PH2.get.cc.EPC
PH2.get.ARLb.EPC
PH2.get.k.EPC
PH2XBAR
PH2XBAR.data
Documented in
PH1.get.cc
PH1XBAR
PH1XBAR.data
PH2.get.ARLb.EPC
PH2.get.cc.EPC
PH2.get.cc.uc
PH2.get.k.EPC
PH2XBAR
PH2XBAR.data
####################################################################################################################################################
# Phase I Xbar Chart Based on Champ and Jones(2004); see also Yao et al. (2017), and Yao and Chakraborti (2018)
####################################################################################################################################################
#use too many functions from package mvtnorm
require(mvtnorm)
#require(adehabitatLT)
#dchi <- adehabitatLT::dchi
#focus on these 3 functions from package mvtnorm
#pmvt <- mvtnorm::pmvt
#qmvt <- mvtnorm::qmvt
#pmvnorm <- mvtnorm::pmvnorm
####################################################################################################################################################
#source('https://raw.githubusercontent.com/bolus123/Statistical-Process-Control/master/MKLswitch.R')
####################################################################################################################################################
W.f <- function(w, n) { # n is sample size
integrand <- function(x, w) {
(pnorm(w + x) - pnorm(x)) ^ (n - 2) * dnorm(w + x) * dnorm(x)
}
n * (n - 1) * integrate(integrand, -Inf, Inf, w = w)$value
}
# moment of W
W.moment <- function(a, n) { # a is order of moment and n is sample size
integrand <- function(w, a, n) {
w ^ a * W.f(w, n)
}
integrand <- Vectorize(integrand, 'w')
integrate(integrand, 0, Inf, a = a, n = n)$value
}
pars.root.finding <- function(k, n, lower = 1e-6) {
root.finding <- function(nu, k, n, ratio) {
ratio / k - nu / 2 / pi * beta(nu / 2, 1 / 2) ^ 2 + 1
}
M <- W.moment(1, n)
M2 <- W.moment(2, n)
V <- M2 - M ^ 2
ratio <- V / M ^ 2
cat('dn2/Vn = ', 1 / ratio, '\n')
nu <- uniroot(root.finding, c(lower, k * n), k = k, n = n, ratio = ratio)$root
lambda <- sqrt(V / k + M ^ 2)
c(nu, lambda)
}
pars.root.finding <- Vectorize(pars.root.finding, c('k', 'n'))
#k <- c(3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 50, 100, 150, 200, 250, 300)
#pars.root.finding(k, 2)
####################################################################################################################################################
#parts of getting the charting constant l
####################################################################################################################################################
d2.f <- function(n) {
integrand <- function(x, n) {
1 - (1 - pnorm(x)) ^ n - pnorm(x) ^ n
}
integrate(integrand, -Inf, Inf, n = n)$value
}
c4.f <- function(nu) sqrt(2 / (nu)) / beta((nu) / 2, 1 / 2) * sqrt(pi) #c4.function
ub.cons.f <- function(nu, ub.option) {
if (ub.option == 'c4') {
ub.cons <- c4.f(nu)
} else if (ub.option == 'd2') {
ub.cons <- d2.f(2)
} else {
ub.cons <- 1
}
ub.cons
}
PH1.corr.f <- function(m, off.diag = - 1 / (m - 1)){
#correlation matrix
crr <- diag(m)
crr[which(crr == 0)] <- off.diag
crr
}
####################################################################################################################################################
#get l using the multivariate t cdf
####################################################################################################################################################
PH1.get.cc.mvt <- function(
m
,nu
,FAP = 0.1
#,Phase1 = TRUE
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,ub.option = TRUE
,maxiter = 10000
){
alternative = '2-sided' #turn off the alternative
corr.P <- PH1.corr.f(m = m, off.diag = off.diag) #get correlation matrix with equal correlations
pu <- 1 - FAP
ub.cons <- ub.cons.f(nu, 'c4')
L <- ifelse(
alternative == '2-sided',
qmvt(pu, df = nu, sigma = corr.P, tail = 'both.tails', maxiter = maxiter)$quantile,
qmvt(pu, df = nu, sigma = corr.P, maxiter = maxiter)$quantile
)
#get L by multivariate T
c.i <- L * ub.cons * sqrt((m - 1) / m) #get c.i
list(l = L, c.i = c.i)
}
####################################################################################################################################################
#get L by multivariate Normal
####################################################################################################################################################
PH1.joint.pdf.mvn.chisq <- function(
Y
,c.i
,m
,nu
,sigma
,lambda = 1
#,alternative = '2-sided'
,ub.cons = 1)
{
alternative = '2-sided' #turn off the alternative
s <- length(Y)
L <- c.i / sqrt((m - 1) / m * nu) * sqrt(Y) / ub.cons * lambda
dpp <- lapply(
1:s,
function(i){
LL <- rep(L[i], m)
ifelse(
alternative == '2-sided',
pmvnorm(lower = -LL, upper = LL, sigma = sigma),
pmvnorm(lower = rep(-Inf, m), upper = LL, sigma = sigma)
)
}
)
dpp <- unlist(dpp)
dpp * dchisq(Y, df = nu)
}
PH1.root.mvn.F <- function(
c.i
, m
, nu
, sigma
, lambda = 1
, pu
#, alternative = '2-sided'
, ub.cons = 1
, subdivisions = 2000
, rel.tol = 1e-2)
{
alternative = '2-sided' #turn off the alternative
pp <- integrate(
PH1.joint.pdf.mvn.chisq,
lower = 0,
upper = Inf,
c.i = c.i,
m = m,
nu = nu,
sigma = sigma,
lambda = lambda,
#alternative = alternative,
ub.cons = ub.cons,
subdivisions = subdivisions,
rel.tol = rel.tol
)$value
pu - pp
}
PH1.get.cc.mvn <- function(
m
,nu
,FAP = 0.1
#,Phase1 = TRUE
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,var.est = 'MSE'
,ub.option = TRUE
,ub.lower = 1e-6
,interval = c(1, 7)
,maxiter = 10000
,subdivisions = 2000
,tol = 1e-2
){
alternative = '2-sided' #turn off the alternative
#The purpose of this function is
#to obtain L and K based on
#MCMC <- FALSE #multivariate normal.
#MCMC part is not available now.
#if (is.null(off.diag)) off.diag <- ifelse(Phase1 == TRUE, - 1 /(m - 1), 1 / (m + 1))
ub.cons <- 1
lambda <- 1
if (var.est == 'MSE') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'c4')
}
} else if (var.est == 'MR') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'd2')
}
nu.lambda <- pars.root.finding(m - 1, 2, lower = ub.lower)
nu <- nu.lambda[1]
lambda <- nu.lambda[2]
} else {
stop('The variance estimation is unknown.')
}
cat('nu:', nu, ', lambda:', lambda, '\n')
corr.P <- PH1.corr.f(m = m, off.diag = off.diag)
pu <- 1 - FAP
c.i <- uniroot(
PH1.root.mvn.F,
interval = interval,
m = m,
nu = nu,
#Y = Y,
sigma = corr.P,
lambda = lambda,
pu = pu,
#alternative = alternative,
ub.cons = ub.cons,
subdivisions = subdivisions,
tol = tol,
rel.tol = tol,
maxiter = maxiter
)$root
L <- c.i / ub.cons * sqrt(m / (m - 1)) * lambda
list(l = L, c.i = c.i)
}
####################################################################################################################################################
#get L by multivariate Normal
####################################################################################################################################################
##############
PH1.rmvn.MC <- function(sim, sigma) {
rmvnorm(sim, sigma = sigma)
}
PH1.joint.pdf.mvn.chisq.MC <- function(
Y
,c.i
,m
,nu
,sigma
,lambda = 1
#,alternative = '2-sided'
,ub.cons = 1
,X
){
alternative = '2-sided' #turn off the alternative
s <- length(Y)
L <- c.i / sqrt((m - 1) / m * nu) * sqrt(Y) / ub.cons * lambda
dpp <- lapply(
1:s,
function(i){
n <- dim(X)[1]
ifelse(
alternative == '2-sided',
mean(rowSums(-L[i] < X & X < L[i]) == m),
mean(rowSums(-Inf < X & X < L[i]) == m)
)
}
)
dpp <- unlist(dpp)
dpp
}
PH1.root.mvn.F.MC <- function(
c.i
, m
, nu
, sigma
, lambda = 1
, pu
#, alternative = '2-sided'
, ub.cons = 1
, X
, Y
){
alternative = '2-sided' #turn off the alternative
pp <- mean(
PH1.joint.pdf.mvn.chisq.MC(
Y = Y,
c.i = c.i,
m = m,
nu = nu,
sigma = sigma,
lambda = lambda,
#alternative = alternative,
ub.cons = ub.cons,
X = X
)
)
cat('ci:', c.i, ' and diff:', pu - pp, '\n')
pu - pp
}
PH1.get.cc.mvn.MC <- function(
m
,nu
,FAP = 0.1
#,Phase1 = TRUE
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,var.est = 'MSE'
,ub.option = TRUE
,ub.lower = 1e-6
,sim.X = 10000
,sim.Y = 10000
,interval = c(1, 7)
,maxiter = 10000
,tol = 1e-2
){
alternative = '2-sided' #turn off the alternative
#The purpose of this function is
#to obtain L and K based on
#MCMC <- FALSE #multivariate normal.
#MCMC part is not available now.
#if (is.null(off.diag)) off.diag <- ifelse(Phase1 == TRUE, - 1 /(m - 1), 1 / (m + 1))
ub.cons <- 1
lambda <- 1
if (var.est == 'MSE') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'c4')
}
} else if (var.est == 'MR') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'd2')
}
nu.lambda <- pars.root.finding(m - 1, 2, lower = ub.lower)
nu <- nu.lambda[1]
lambda <- nu.lambda[2]
} else {
stop('The variance estimation is unknown.')
}
cat('nu:', nu, ', lambda:', lambda, '\n')
corr.P <- PH1.corr.f(m = m, off.diag = off.diag)
X <- PH1.rmvn.MC(sim = sim.X, sigma = corr.P)
Y <- rchisq(sim.Y, nu)
pu <- 1 - FAP
c.i <- uniroot(
PH1.root.mvn.F.MC,
interval = interval,
m = m,
nu = nu,
#Y = Y,
sigma = corr.P,
lambda = lambda,
pu = pu,
#alternative = alternative,
ub.cons = ub.cons,
X = X,
Y = Y,
tol = tol,
maxiter = maxiter
)$root
L <- c.i / ub.cons * sqrt(m / (m - 1)) * lambda
list(l = L, c.i = c.i)
}
####################################################################################################################################################
PH1.get.cc <- function(
m
,nu = m - 1
,FAP = 0.1
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,ub.option = TRUE
,maxiter = 10000
,method = 'direct'
,var.est = 'MSE'
,ub.lower = 1e-6
,indirect.interval = c(1, 7)
,indirect.subdivisions = 100L
,indirect.tol = .Machine$double.eps^0.25
){
alternative = '2-sided' #turn off the alternative
#The purpose of this function is to obtain L and K
#by multivariate T or multivariate normal.
#Multivariate normal is so time-consuming
#that I do not recommend.
is.int <- ifelse(nu == round(nu), 1, 0)
if (method == 'direct') {
if (var.est == 'MSE') {
if (is.int == 1) {
PH1.get.cc.mvt(
m = m
,nu = nu
,FAP = FAP
#,Phase1 = Phase1
,off.diag = off.diag
#,alternative = alternative
,ub.option = ub.option
,maxiter = maxiter
)
} else if (is.int == 0) {
stop('Nu is not an integer. Please use the indirect method instead.')
}
} else {
stop('The variance estimation must be MSE for the direct method. Please use the indirect method instead.')
}
} else if (method == 'indirect') {
if (is.int == 1 & var.est == 'MSE') {
cat('Nu is an integer. Using the indirect method may slow the computation process down.', '\n')
}
PH1.get.cc.mvn(
m = m
,nu = nu
,FAP = FAP
#,Phase1 = Phase1
,off.diag = off.diag
#,alternative = alternative
,var.est = var.est
,ub.option = ub.option
,ub.lower = ub.lower
,interval = indirect.interval
#,maxsim = indirect.maxsim
,subdivisions = indirect.subdivisions
,maxiter = maxiter
,tol = indirect.tol
)
} else {
stop('Unknown method. Please select the direct method or the indirect method.')
}
}
####################################################################################################################################################
#reverse the process
####################################################################################################################################################
#only support the direct method
#PH1.get.FAP0 <- function(
# c.i
# ,k
# ,nu
# ,off.diag = -1/(k - 1)
# ,alternative = '2-sided'
# ,c4.option = TRUE
# ,maxiter = 10000
# ,indirect.subdivisions = 100L
# ,indirect.tol = .Machine$double.eps^0.25
#
#){
#
# #method <- 'direct'
#
# #Phase1 <- TRUE
#
# is.int <- ifelse(nu == round(nu), 1, 0)
#
# #if (is.null(off.diag)) off.diag <- ifelse(Phase1 == TRUE, - 1 /(m - 1), 1 / (m + 1))
#
# corr.P <- PH1.corr.f(k = k, off.diag = off.diag)
#
# corr.par <- corr.par.f(nu, c4.option)
#
# L <- c.i / corr.par * sqrt(k / (k - 1))
#
#
# if (method == 'direct' & is.int == 1) { #using multivariate T to obtain L and K
#
# ll <- rep(L, k)
#
# ifelse(
# alternative == '2-sided',
# 1 - pmvt(lower = -ll, upper = ll, df = nu, corr = corr.P, algorithm = TVPACK, abseps = 1e-12),
# 1 - pmvt(lower = -Inf, upper = ll, df = nu, corr = corr.P, algorithm = TVPACK, abseps = 1e-12)
# )
#
# }
#
#
#
#
#}
####################################################################################################################################################
#Build Control Chart
####################################################################################################################################################
PH1XBAR <- function(
X
,c.i = NULL
,FAP = 0.1
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,model = 'ANOVA-based'
,ub.option = 'c4'
,plot.option = TRUE
,maxiter = 10000
,method = 'direct'
,indirect.interval = c(1, 7)
,indirect.subdivisions = 100L
,indirect.tol = .Machine$double.eps^0.25
) {
alternative = '2-sided' #turn off the alternative
m <- dim(X)[1]
n <- dim(X)[2]
X.bar <- rowMeans(X)
X.bar.bar <- mean(X)
if (model == 'ANOVA-based') {
nu <- m - 1
ub.cons <- ub.cons.f(nu, ub.option)
sigma.v <- sqrt(var(X.bar)) / ub.cons
} else if (model == 'basic') {
nu <- m * (n - 1)
ub.cons <- ub.cons.f(nu, ub.option)
sigma.v <- sqrt(sum(apply(X, 1, var)) / m) / ub.cons / sqrt(n)
} else {
stop("need to specify whether it is based on the ANOVA-based model or others")
}
if (is.null(c.i)) {
c.i <- PH1.get.cc(
m = m
,nu = nu
,FAP = FAP
,off.diag = off.diag
#,alternative = alternative
,ub.option = ub.option
,maxiter = maxiter
,method = method
,indirect.interval = indirect.interval
,indirect.subdivisions = indirect.subdivisions
,indirect.tol = indirect.tol
)$c.i
} else {
c.i <- c.i
}
LCL <- X.bar.bar - c.i * sigma.v
UCL <- X.bar.bar + c.i * sigma.v
if (plot.option == TRUE){
plot(c(1, m), c(LCL, UCL), xaxt = "n", xlab = 'Subgroup', ylab = 'Sample Mean', type = 'n')
axis(side = 1, at = 1:m)
points(1:m, X.bar, type = 'o')
points(c(-1, m + 2), c(LCL, LCL), type = 'l', col = 'red')
points(c(-1, m + 2), c(UCL, UCL), type = 'l', col = 'red')
points(c(-1, m + 2), c(X.bar.bar, X.bar.bar), type = 'l', col = 'blue')
text(round(m * 0.8), UCL, paste('UCL = ', round(UCL, 4)), pos = 1)
text(round(m * 0.8), LCL, paste('LCL = ', round(LCL, 4)), pos = 3)
}
list(CL = X.bar.bar, sigma = sigma.v, PH1.cc = c.i, m = m, nu = nu, LCL = LCL, UCL = UCL, CS = X.bar)
}
####################################################################################################################################################
#Example about table 6.3 (Montgomery, 2013)
####################################################################################################################################################
PH1XBAR.data <- function(){
matrix(c(
74.03 , 74.002 , 74.019 , 73.992 , 74.008 ,
73.995 , 73.992 , 74.001 , 74.011 , 74.004 ,
73.988 , 74.024 , 74.021 , 74.005 , 74.002 ,
74.002 , 73.996 , 73.993 , 74.015 , 74.009 ,
73.992 , 74.007 , 74.015 , 73.989 , 74.014 ,
74.009 , 73.994 , 73.997 , 73.985 , 73.993 ,
73.995 , 74.006 , 73.994 , 74 , 74.005 ,
73.985 , 74.003 , 73.993 , 74.015 , 73.988 ,
74.008 , 73.995 , 74.009 , 74.005 , 74.004 ,
73.998 , 74 , 73.99 , 74.007 , 73.995 ,
73.994 , 73.998 , 73.994 , 73.995 , 73.99 ,
74.004 , 74 , 74.007 , 74 , 73.996 ,
73.983 , 74.002 , 73.998 , 73.997 , 74.012 ,
74.006 , 73.967 , 73.994 , 74 , 73.984 ,
74.012 , 74.014 , 73.998 , 73.999 , 74.007 ,
74 , 73.984 , 74.005 , 73.998 , 73.996 ,
73.994 , 74.012 , 73.986 , 74.005 , 74.007 ,
74.006 , 74.01 , 74.018 , 74.003 , 74 ,
73.984 , 74.002 , 74.003 , 74.005 , 73.997 ,
74 , 74.01 , 74.013 , 74.02 , 74.003 ,
73.982 , 74.001 , 74.015 , 74.005 , 73.996 ,
74.004 , 73.999 , 73.99 , 74.006 , 74.009 ,
74.01 , 73.989 , 73.99 , 74.009 , 74.014 ,
74.015 , 74.008 , 73.993 , 74 , 74.01 ,
73.982 , 73.984 , 73.995 , 74.017 , 74.013
), ncol = 5, byrow = T)
}
####################################################################################################################################################
# Phase II Xbar chart. Please see Yao and Chakraborti (2018)
####################################################################################################################################################
PH2.inner.normal <- function(Z, Y, c.ii, k, nu = k - 1, c4.option = TRUE){
corr.par <- corr.par.f(nu, c4.option)
qn <- Z / sqrt(k) + c.ii / corr.par / sqrt(nu) * sqrt(Y)
pnorm(qn)
}
PH2.CFAR.intgrand <- function(u, v, c.ii, k, nu = k - 1, c4.option = TRUE){
Z <- qnorm(u)
Y <- qchisq(v, nu)
p1 <- PH2.inner.normal(Z, Y, c.ii, k, nu, c4.option)
p2 <- PH2.inner.normal(Z, Y, -c.ii, k, nu, c4.option)
1 - p1 + p2
}
PH2.CARL.intgrand <- function(u, v, c.ii, k, nu = k - 1, c4.option = TRUE) 1 / PH2.CFAR.intgrand(u, v, c.ii, k, nu, c4.option)
####################################################################################################################################################
PH2.get.cc.uc <- function(ARL0, k, nu = k - 1, c4.option = TRUE, interval = c(1, 10), u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
PH2.root.finding.uc <- function(c.ii, k, nu = k - 1, c4.option = TRUE, ARL0, u, v) {
ARL0 - mean(PH2.CARL.intgrand(u, v, c.ii, k, nu, c4.option))
}
k <- k
rt <- uniroot(PH2.root.finding.uc, interval = interval, k = k, nu = nu, c4.option = c4.option, ARL0 = ARL0, u = u, v = v, tol = tol)$root
rt
}
####################################################################################################################################################
PH2.root.finding.EPC <- function(p, k, nu = k - 1, c.ii, ARLb, c4.option = TRUE, u = runif(100000), v = runif(100000)){
1 - p - mean(PH2.CARL.intgrand(u, v, c.ii, k, nu, c4.option) >= ARLb)
}
PH2.get.cc.EPC <- function(p, k, nu = k - 1, eps = 0.1, ARL0 = 370, c4.option = TRUE, interval = c(1, 10), u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
ARLb <- (1 - eps) * ARL0
rt <- uniroot(
PH2.root.finding.EPC,
interval = interval,
p = p,
k = k,
nu = nu,
ARLb = ARLb,
c4.option = c4.option,
u = u,
v = v,
tol = tol
)$root
rt
}
#rnum <- 100000
#u <- runif(rnum)
#v <- runif(rnum)
#debug(PH2.get.c.ii.EPC)
#PH2.get.c.ii.EPC(p = 0.05, k = 100, eps = 0.1, ARL0 = 500, c4.option = TRUE, interval = c(1, 5), u = u, v = v)
PH2.get.ARLb.EPC <- function(p, k, nu = k - 1, c.ii, ARL0 = NULL, c4.option = TRUE, interval = c(1, 1000), u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
rt <- uniroot(
PH2.root.finding.EPC,
interval = interval,
p = p,
k = k,
nu = nu,
c.ii = c.ii,
c4.option = c4.option,
u = u,
v = v,
tol = tol
)$root
if (is.null(ARL0)) {
list(ARLb = rt, eps = NULL, ARL0 = NULL)
} else {
list(ARLb = rt, eps = 1 - rt / ARL0, ARL0 = ARL0)
}
}
#PH2.get.ARLb.EPC(p = 0.1, k = 100, c.ii = 3, c4.option = TRUE, u = u, v = v)
#PH2.get.ARLb.EPC(p = 0.1, k = 100, c.ii = 3, ARL0 = 370, c4.option = TRUE, u = u, v = v)
PH2.get.k.EPC <- function(p, c.ii, eps = 0.1, ARL0 = 370, c4.option = TRUE, interval = c(1000, 3000), model = 'ANOVA-based', n = 10, u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
if (model == 'ANOVA-based') {
PH2.root.finding.k.EPC <- function(p, k, c.ii, ARLb, c4.option = TRUE, n = 10, u, v){
nu <- k - 1
PH2.root.finding.EPC(p, k, nu, c.ii, ARLb, c4.option, u, v)
}
} else if (model == 'basic') {
PH2.root.finding.k.EPC <- function(p, k, c.ii, ARLb, c4.option = TRUE, n = 10, u, v){
nu <- (n - 1) * k
PH2.root.finding.EPC(p, k, nu, c.ii, ARLb, c4.option, u, v)
}
} else {
stop("need to specify whether it is based on the ANOVA-based model or others")
}
ARLb <- (1 - eps) * ARL0
rt <- uniroot(
PH2.root.finding.k.EPC,
interval = interval,
p = p,
c.ii = c.ii,
ARLb = ARLb,
c4.option = c4.option,
n = n,
u = u,
v = v,
tol = tol
)$root
list(k = rt, PH1.sample.size = rt * n)
}
####################################################################################################################################################
PH2XBAR <- function(
X
,PH1.info = list(X = NULL, mu = NULL, sigma = NULL, k = NULL, n = NULL, model = "ANOVA-based")
,c.ii.info = list(c.ii = NULL, method = 'UC', ARL0 = 370, p = 0.05, eps = 0.1, interval.c.ii.UC = c(1, 3.2), interval.c.ii.EPC = c(1, 10), UC.tol = .Machine$double.eps^0.25, EPC.tol = .Machine$double.eps^0.25)
#,alternative = '2-sided'
,c4.option = TRUE
,plot.option = TRUE
,maxsim = 100000
) {
alternative = '2-sided' #turn off the alternative
k.ii <- dim(X)[1]
X.bar <- rowMeans(X)
if (is.null(PH1.info$X)){
k <- PH1.info$k
if (PH1.info$model == 'ANOVA-based') {
nu <- k - 1
} else if (PH1.info$model == 'basic') {
if (is.null(PH1.info$n)) {
stop('Subgroup size needs to be defined')
} else {
n <- PH1.info$n
nu <- k * (n - 1)
}
} else {
stop('The model in Phase 1 needs to be defined')
}
corr.par <- corr.par.f(nu, c4.option)
X.bar.bar <- PH1.info$mu
sigma.v <- PH1.info$sigma / corr.par
} else {
PH1.chart <- PH1XBAR(X = PH1.info$X, c.i = 1, model = PH1.info$model, plot.option = FALSE)
k <- PH1.chart$k
nu <- PH1.chart$nu
X.bar.bar <- PH1.chart$CL
sigma.v <- PH1.chart$sigma
}
u <- runif(maxsim)
v <- runif(maxsim)
UC <- NULL
EPC <- NULL
n.c.ii <- 0
UC.flag <- 0
EPC.flag <- 0
if (is.null(c.ii.info$c.ii)) {
if (c.ii.info$method == 'both') {
UC.flag <- 1
EPC.flag <- 1
} else if (c.ii.info$method == 'UC') {
UC.flag <- 1
} else if (c.ii.info$method == 'EPC'){
EPC.flag <- 1
} else {
stop("Please check the method of obtaining charting constants in Phase 2")
}
if (UC.flag == 1) {
n.c.ii <- n.c.ii + 1
UC <- PH2.get.cc.uc(
ARL0 = c.ii.info$ARL0
, k = k
, nu = nu
, c4.option = c4.option
, interval = c.ii.info$interval.c.ii.UC
, u = u
, v = v
, tol = c.ii.info$UC.tol
)
}
if (EPC.flag == 1){
n.c.ii <- n.c.ii + 1
EPC <- PH2.get.cc.EPC(
p = c.ii.info$p
, k = k
, nu = nu
, eps = c.ii.info$eps
, ARL0 = c.ii.info$ARL0
, c4.option = c4.option
, interval = c.ii.info$interval.c.ii.EPC
, u = u
, v = v
, tol = c.ii.info$EPC.tol
)
}
c.ii <- list(UC = UC, EPC = EPC)
} else {
n.c.ii <- length(c.ii.info$c.ii)
c.ii <- c.ii.info$c.ii
}
LCL <- X.bar.bar - unlist(c.ii) * sigma.v
UCL <- X.bar.bar + unlist(c.ii) * sigma.v
if (plot.option == TRUE){
plot(c(1, k.ii), c(min(c(LCL, X.bar)), max(c(UCL, X.bar))), xaxt = "n", xlab = 'Subgroup', ylab = 'Sample Mean', type = 'n')
points(1:k.ii, X.bar, type = 'o')
axis(side = 1, at = 1:k.ii)
points(c(-1, k.ii + 2), c(X.bar.bar, X.bar.bar), type = 'l', col = 'blue')
for (ii in 1:n.c.ii){
points(c(-1, k.ii + 2), c(LCL[ii], LCL[ii]), type = 'l', col = 'red')
points(c(-1, k.ii + 2), c(UCL[ii], UCL[ii]), type = 'l', col = 'red')
if (is.null(c.ii.info$c.ii)) {
lab.in <- paste('(', names(c.ii)[ii], ')', sep = '')
} else {
lab.in <- NULL
}
ulab <- paste('UCL', lab.in, ' = ', round(UCL[ii], 4), sep = '')
llab <- paste('LCL', lab.in, ' = ', round(LCL[ii], 4), sep = '')
text(round(k.ii * 0.8), UCL[ii], ulab, pos = 1)
text(round(k.ii * 0.8), LCL[ii], llab, pos = 3)
}
}
list(CL = X.bar.bar, sigma = sigma.v, k = k, nu = nu, PH2.cc = unlist(c.ii), LCL = LCL, UCL = UCL, CS = X.bar)
}
####################################################################################################################################################
PH2XBAR.data <- function(){
matrix(
c(
240, 243, 250, 253, 248,
238, 242, 245, 251, 247,
239, 242, 246, 250, 248,
235, 237, 246, 249, 246,
240, 241, 246, 247, 249,
240, 243, 244, 248, 245,
240, 243, 244, 249, 246,
245, 250, 250, 247, 248,
238, 240, 245, 248, 246,
240, 242, 246, 249, 248,
240, 243, 246, 250, 248,
241, 245, 243, 247, 245,
247, 245, 255, 250, 249,
237, 239, 243, 247, 246,
242, 244, 245, 248, 245,
237, 239, 242, 247, 245,
242, 244, 246, 251, 248,
243, 245, 247, 252, 249,
243, 245, 248, 251, 250,
244, 246, 246, 250, 246,
241, 239, 244, 250, 246,
242, 245, 248, 251, 249,
242, 245, 248, 243, 246,
241, 244, 245, 249, 247,
236, 239, 241, 246, 242,
243, 246, 247, 252, 247,
241, 243, 245, 248, 246,
239, 240, 242, 243, 244,
239, 240, 250, 252, 250,
241, 243, 249, 255, 253
)
,nrow = 30
,ncol = 5
,byrow = TRUE
)
}
#PH2.get.k.EPC(p = 0.05, c.ii = 3, eps = 0.2, ARL0 = 370, c4.option = TRUE, interval = c(100, 400), type = 'ANOVA-based', subgroup.size = 10, u = u, v = v)
#
#m <- dim(x)[1]
#n <- dim(x)[2]
#
#x.bar.bar <- mean(x)
#
#sigma.v <- sqrt(sum(apply(x, 1, var)) / m) / c4.f(m * (n - 1))
#
#k <- get.cc(m, m * (n - 1), 0.05)$K
#
#LCL <- x.bar.bar - k * sigma.v / sqrt(n)
#UCL <- x.bar.bar + k * sigma.v / sqrt(n)
#
#
#x.bar <- rowMeans(x)
#
#plot(c(1, m), c(LCL * 49997/50000, UCL * 50003/50000), xaxt = "n", xlab = 'Subgroup', ylab = 'Sample Means', type = 'n')
#
#axis(side = 1, at = 1:25)
#
#points(1:m, x.bar, type = 'o')
#points(c(-1, m + 2), c(LCL, LCL), type = 'l', col = 'red')
#points(c(-1, m + 2), c(UCL, UCL), type = 'l', col = 'red')
#points(c(-1, m + 2), c(x.bar.bar, x.bar.bar), type = 'l', col = 'blue')
#text(20, UCL * 50001/50000, paste('UCL = ', round(UCL, 4)))
#text(20, LCL * 49999/50000, paste('LCL = ', round(LCL, 4)))
#
#pos <- c(rep(c(3, 1), 5), 1, 3, 3, 1, 3, 1, 1, 3, 1, 3, 1, 3, 1, 3, 1) #rep(c(1, 3), 8))[-26]
#
#text(1:m, x.bar, round(x.bar, 4), cex = 0.7, pos = pos)
#
#
####################################################################################################################################################
#Example about how to get K and L by FAP0
####################################################################################################################################################
#get.cc(
# 20
# ,80
# ,FAP = 0.1
# ,off.diag = NULL
# ,alternative = '2-sided'
# ,maxiter = 10000
# ,method = 'direct'
#)
#get.cc(
# 20
# ,80
# ,FAP = 0.1
# ,off.diag = NULL
# ,alternative = '2-sided'
# ,maxiter = 10000
# ,method = 'indirect'
#)
####################################################################################################################################################
#Example about how to get FAP by K
####################################################################################################################################################
#
#k <- c(
#22.84794
#,10.52163
#,7.386245
#,6.057107
#,5.340985
#,4.898076
#,4.598652
#,4.383316
#,3.842428
#,3.619871
#,3.498807
#,3.422752
#,3.280788
#,3.182391
#,3.151007
#,3.135567
#,3.126384
#,3.120294
#
#
#
#
#
#)
#
#m.seq <- c(
# 3
# ,4
# ,5
# ,6
# ,7
# ,8
# ,9
# ,10
# ,15
# ,20
# ,25
# ,30
# ,50
# ,100
# ,150
# ,200
# ,250
# ,300
#
#)
#
#
#i <- 0
#
#record <- rep(NA, 18)
#
#for (m in m.seq){
#
# i <- i + 1
#
# nu <- m - 1
#
# record[i] <- get.FAP0(k[i], m, nu)
#
#}
#
#as.matrix(record)
bolus123/PH1AND2XBAR documentation built on May 7, 2019, 10:56 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions W.f W.moment pars.root.finding d2.f c4.f ub.cons.f PH1.corr.f PH1.get.cc.mvt PH1.joint.pdf.mvn.chisq PH1.root.mvn.F PH1.get.cc.mvn PH1.rmvn.MC PH1.joint.pdf.mvn.chisq.MC PH1.root.mvn.F.MC PH1.get.cc.mvn.MC PH1.get.cc PH1XBAR PH1XBAR.data PH2.inner.normal PH2.CFAR.intgrand PH2.CARL.intgrand PH2.get.cc.uc PH2.root.finding.EPC PH2.get.cc.EPC PH2.get.ARLb.EPC PH2.get.k.EPC PH2XBAR PH2XBAR.data
Documented in PH1.get.cc PH1XBAR PH1XBAR.data PH2.get.ARLb.EPC PH2.get.cc.EPC PH2.get.cc.uc PH2.get.k.EPC PH2XBAR PH2XBAR.data
####################################################################################################################################################
# Phase I Xbar Chart Based on Champ and Jones(2004); see also Yao et al. (2017), and Yao and Chakraborti (2018)
####################################################################################################################################################
#use too many functions from package mvtnorm
require(mvtnorm)
#require(adehabitatLT)
#dchi <- adehabitatLT::dchi
#focus on these 3 functions from package mvtnorm
#pmvt <- mvtnorm::pmvt
#qmvt <- mvtnorm::qmvt
#pmvnorm <- mvtnorm::pmvnorm
####################################################################################################################################################
#source('https://raw.githubusercontent.com/bolus123/Statistical-Process-Control/master/MKLswitch.R')
####################################################################################################################################################
W.f <- function(w, n) { # n is sample size
integrand <- function(x, w) {
(pnorm(w + x) - pnorm(x)) ^ (n - 2) * dnorm(w + x) * dnorm(x)
}
n * (n - 1) * integrate(integrand, -Inf, Inf, w = w)$value
}
# moment of W
W.moment <- function(a, n) { # a is order of moment and n is sample size
integrand <- function(w, a, n) {
w ^ a * W.f(w, n)
}
integrand <- Vectorize(integrand, 'w')
integrate(integrand, 0, Inf, a = a, n = n)$value
}
pars.root.finding <- function(k, n, lower = 1e-6) {
root.finding <- function(nu, k, n, ratio) {
ratio / k - nu / 2 / pi * beta(nu / 2, 1 / 2) ^ 2 + 1
}
M <- W.moment(1, n)
M2 <- W.moment(2, n)
V <- M2 - M ^ 2
ratio <- V / M ^ 2
cat('dn2/Vn = ', 1 / ratio, '\n')
nu <- uniroot(root.finding, c(lower, k * n), k = k, n = n, ratio = ratio)$root
lambda <- sqrt(V / k + M ^ 2)
c(nu, lambda)
}
pars.root.finding <- Vectorize(pars.root.finding, c('k', 'n'))
#k <- c(3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 50, 100, 150, 200, 250, 300)
#pars.root.finding(k, 2)
####################################################################################################################################################
#parts of getting the charting constant l
####################################################################################################################################################
d2.f <- function(n) {
integrand <- function(x, n) {
1 - (1 - pnorm(x)) ^ n - pnorm(x) ^ n
}
integrate(integrand, -Inf, Inf, n = n)$value
}
c4.f <- function(nu) sqrt(2 / (nu)) / beta((nu) / 2, 1 / 2) * sqrt(pi) #c4.function
ub.cons.f <- function(nu, ub.option) {
if (ub.option == 'c4') {
ub.cons <- c4.f(nu)
} else if (ub.option == 'd2') {
ub.cons <- d2.f(2)
} else {
ub.cons <- 1
}
ub.cons
}
PH1.corr.f <- function(m, off.diag = - 1 / (m - 1)){
#correlation matrix
crr <- diag(m)
crr[which(crr == 0)] <- off.diag
crr
}
####################################################################################################################################################
#get l using the multivariate t cdf
####################################################################################################################################################
PH1.get.cc.mvt <- function(
m
,nu
,FAP = 0.1
#,Phase1 = TRUE
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,ub.option = TRUE
,maxiter = 10000
){
alternative = '2-sided' #turn off the alternative
corr.P <- PH1.corr.f(m = m, off.diag = off.diag) #get correlation matrix with equal correlations
pu <- 1 - FAP
ub.cons <- ub.cons.f(nu, 'c4')
L <- ifelse(
alternative == '2-sided',
qmvt(pu, df = nu, sigma = corr.P, tail = 'both.tails', maxiter = maxiter)$quantile,
qmvt(pu, df = nu, sigma = corr.P, maxiter = maxiter)$quantile
)
#get L by multivariate T
c.i <- L * ub.cons * sqrt((m - 1) / m) #get c.i
list(l = L, c.i = c.i)
}
####################################################################################################################################################
#get L by multivariate Normal
####################################################################################################################################################
PH1.joint.pdf.mvn.chisq <- function(
Y
,c.i
,m
,nu
,sigma
,lambda = 1
#,alternative = '2-sided'
,ub.cons = 1)
{
alternative = '2-sided' #turn off the alternative
s <- length(Y)
L <- c.i / sqrt((m - 1) / m * nu) * sqrt(Y) / ub.cons * lambda
dpp <- lapply(
1:s,
function(i){
LL <- rep(L[i], m)
ifelse(
alternative == '2-sided',
pmvnorm(lower = -LL, upper = LL, sigma = sigma),
pmvnorm(lower = rep(-Inf, m), upper = LL, sigma = sigma)
)
}
)
dpp <- unlist(dpp)
dpp * dchisq(Y, df = nu)
}
PH1.root.mvn.F <- function(
c.i
, m
, nu
, sigma
, lambda = 1
, pu
#, alternative = '2-sided'
, ub.cons = 1
, subdivisions = 2000
, rel.tol = 1e-2)
{
alternative = '2-sided' #turn off the alternative
pp <- integrate(
PH1.joint.pdf.mvn.chisq,
lower = 0,
upper = Inf,
c.i = c.i,
m = m,
nu = nu,
sigma = sigma,
lambda = lambda,
#alternative = alternative,
ub.cons = ub.cons,
subdivisions = subdivisions,
rel.tol = rel.tol
)$value
pu - pp
}
PH1.get.cc.mvn <- function(
m
,nu
,FAP = 0.1
#,Phase1 = TRUE
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,var.est = 'MSE'
,ub.option = TRUE
,ub.lower = 1e-6
,interval = c(1, 7)
,maxiter = 10000
,subdivisions = 2000
,tol = 1e-2
){
alternative = '2-sided' #turn off the alternative
#The purpose of this function is
#to obtain L and K based on
#MCMC <- FALSE #multivariate normal.
#MCMC part is not available now.
#if (is.null(off.diag)) off.diag <- ifelse(Phase1 == TRUE, - 1 /(m - 1), 1 / (m + 1))
ub.cons <- 1
lambda <- 1
if (var.est == 'MSE') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'c4')
}
} else if (var.est == 'MR') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'd2')
}
nu.lambda <- pars.root.finding(m - 1, 2, lower = ub.lower)
nu <- nu.lambda[1]
lambda <- nu.lambda[2]
} else {
stop('The variance estimation is unknown.')
}
cat('nu:', nu, ', lambda:', lambda, '\n')
corr.P <- PH1.corr.f(m = m, off.diag = off.diag)
pu <- 1 - FAP
c.i <- uniroot(
PH1.root.mvn.F,
interval = interval,
m = m,
nu = nu,
#Y = Y,
sigma = corr.P,
lambda = lambda,
pu = pu,
#alternative = alternative,
ub.cons = ub.cons,
subdivisions = subdivisions,
tol = tol,
rel.tol = tol,
maxiter = maxiter
)$root
L <- c.i / ub.cons * sqrt(m / (m - 1)) * lambda
list(l = L, c.i = c.i)
}
####################################################################################################################################################
#get L by multivariate Normal
####################################################################################################################################################
##############
PH1.rmvn.MC <- function(sim, sigma) {
rmvnorm(sim, sigma = sigma)
}
PH1.joint.pdf.mvn.chisq.MC <- function(
Y
,c.i
,m
,nu
,sigma
,lambda = 1
#,alternative = '2-sided'
,ub.cons = 1
,X
){
alternative = '2-sided' #turn off the alternative
s <- length(Y)
L <- c.i / sqrt((m - 1) / m * nu) * sqrt(Y) / ub.cons * lambda
dpp <- lapply(
1:s,
function(i){
n <- dim(X)[1]
ifelse(
alternative == '2-sided',
mean(rowSums(-L[i] < X & X < L[i]) == m),
mean(rowSums(-Inf < X & X < L[i]) == m)
)
}
)
dpp <- unlist(dpp)
dpp
}
PH1.root.mvn.F.MC <- function(
c.i
, m
, nu
, sigma
, lambda = 1
, pu
#, alternative = '2-sided'
, ub.cons = 1
, X
, Y
){
alternative = '2-sided' #turn off the alternative
pp <- mean(
PH1.joint.pdf.mvn.chisq.MC(
Y = Y,
c.i = c.i,
m = m,
nu = nu,
sigma = sigma,
lambda = lambda,
#alternative = alternative,
ub.cons = ub.cons,
X = X
)
)
cat('ci:', c.i, ' and diff:', pu - pp, '\n')
pu - pp
}
PH1.get.cc.mvn.MC <- function(
m
,nu
,FAP = 0.1
#,Phase1 = TRUE
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,var.est = 'MSE'
,ub.option = TRUE
,ub.lower = 1e-6
,sim.X = 10000
,sim.Y = 10000
,interval = c(1, 7)
,maxiter = 10000
,tol = 1e-2
){
alternative = '2-sided' #turn off the alternative
#The purpose of this function is
#to obtain L and K based on
#MCMC <- FALSE #multivariate normal.
#MCMC part is not available now.
#if (is.null(off.diag)) off.diag <- ifelse(Phase1 == TRUE, - 1 /(m - 1), 1 / (m + 1))
ub.cons <- 1
lambda <- 1
if (var.est == 'MSE') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'c4')
}
} else if (var.est == 'MR') {
if (ub.option == TRUE) {
ub.cons <- ub.cons.f(nu, 'd2')
}
nu.lambda <- pars.root.finding(m - 1, 2, lower = ub.lower)
nu <- nu.lambda[1]
lambda <- nu.lambda[2]
} else {
stop('The variance estimation is unknown.')
}
cat('nu:', nu, ', lambda:', lambda, '\n')
corr.P <- PH1.corr.f(m = m, off.diag = off.diag)
X <- PH1.rmvn.MC(sim = sim.X, sigma = corr.P)
Y <- rchisq(sim.Y, nu)
pu <- 1 - FAP
c.i <- uniroot(
PH1.root.mvn.F.MC,
interval = interval,
m = m,
nu = nu,
#Y = Y,
sigma = corr.P,
lambda = lambda,
pu = pu,
#alternative = alternative,
ub.cons = ub.cons,
X = X,
Y = Y,
tol = tol,
maxiter = maxiter
)$root
L <- c.i / ub.cons * sqrt(m / (m - 1)) * lambda
list(l = L, c.i = c.i)
}
####################################################################################################################################################
PH1.get.cc <- function(
m
,nu = m - 1
,FAP = 0.1
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,ub.option = TRUE
,maxiter = 10000
,method = 'direct'
,var.est = 'MSE'
,ub.lower = 1e-6
,indirect.interval = c(1, 7)
,indirect.subdivisions = 100L
,indirect.tol = .Machine$double.eps^0.25
){
alternative = '2-sided' #turn off the alternative
#The purpose of this function is to obtain L and K
#by multivariate T or multivariate normal.
#Multivariate normal is so time-consuming
#that I do not recommend.
is.int <- ifelse(nu == round(nu), 1, 0)
if (method == 'direct') {
if (var.est == 'MSE') {
if (is.int == 1) {
PH1.get.cc.mvt(
m = m
,nu = nu
,FAP = FAP
#,Phase1 = Phase1
,off.diag = off.diag
#,alternative = alternative
,ub.option = ub.option
,maxiter = maxiter
)
} else if (is.int == 0) {
stop('Nu is not an integer. Please use the indirect method instead.')
}
} else {
stop('The variance estimation must be MSE for the direct method. Please use the indirect method instead.')
}
} else if (method == 'indirect') {
if (is.int == 1 & var.est == 'MSE') {
cat('Nu is an integer. Using the indirect method may slow the computation process down.', '\n')
}
PH1.get.cc.mvn(
m = m
,nu = nu
,FAP = FAP
#,Phase1 = Phase1
,off.diag = off.diag
#,alternative = alternative
,var.est = var.est
,ub.option = ub.option
,ub.lower = ub.lower
,interval = indirect.interval
#,maxsim = indirect.maxsim
,subdivisions = indirect.subdivisions
,maxiter = maxiter
,tol = indirect.tol
)
} else {
stop('Unknown method. Please select the direct method or the indirect method.')
}
}
####################################################################################################################################################
#reverse the process
####################################################################################################################################################
#only support the direct method
#PH1.get.FAP0 <- function(
# c.i
# ,k
# ,nu
# ,off.diag = -1/(k - 1)
# ,alternative = '2-sided'
# ,c4.option = TRUE
# ,maxiter = 10000
# ,indirect.subdivisions = 100L
# ,indirect.tol = .Machine$double.eps^0.25
#
#){
#
# #method <- 'direct'
#
# #Phase1 <- TRUE
#
# is.int <- ifelse(nu == round(nu), 1, 0)
#
# #if (is.null(off.diag)) off.diag <- ifelse(Phase1 == TRUE, - 1 /(m - 1), 1 / (m + 1))
#
# corr.P <- PH1.corr.f(k = k, off.diag = off.diag)
#
# corr.par <- corr.par.f(nu, c4.option)
#
# L <- c.i / corr.par * sqrt(k / (k - 1))
#
#
# if (method == 'direct' & is.int == 1) { #using multivariate T to obtain L and K
#
# ll <- rep(L, k)
#
# ifelse(
# alternative == '2-sided',
# 1 - pmvt(lower = -ll, upper = ll, df = nu, corr = corr.P, algorithm = TVPACK, abseps = 1e-12),
# 1 - pmvt(lower = -Inf, upper = ll, df = nu, corr = corr.P, algorithm = TVPACK, abseps = 1e-12)
# )
#
# }
#
#
#
#
#}
####################################################################################################################################################
#Build Control Chart
####################################################################################################################################################
PH1XBAR <- function(
X
,c.i = NULL
,FAP = 0.1
,off.diag = -1/(m - 1)
#,alternative = '2-sided'
,model = 'ANOVA-based'
,ub.option = 'c4'
,plot.option = TRUE
,maxiter = 10000
,method = 'direct'
,indirect.interval = c(1, 7)
,indirect.subdivisions = 100L
,indirect.tol = .Machine$double.eps^0.25
) {
alternative = '2-sided' #turn off the alternative
m <- dim(X)[1]
n <- dim(X)[2]
X.bar <- rowMeans(X)
X.bar.bar <- mean(X)
if (model == 'ANOVA-based') {
nu <- m - 1
ub.cons <- ub.cons.f(nu, ub.option)
sigma.v <- sqrt(var(X.bar)) / ub.cons
} else if (model == 'basic') {
nu <- m * (n - 1)
ub.cons <- ub.cons.f(nu, ub.option)
sigma.v <- sqrt(sum(apply(X, 1, var)) / m) / ub.cons / sqrt(n)
} else {
stop("need to specify whether it is based on the ANOVA-based model or others")
}
if (is.null(c.i)) {
c.i <- PH1.get.cc(
m = m
,nu = nu
,FAP = FAP
,off.diag = off.diag
#,alternative = alternative
,ub.option = ub.option
,maxiter = maxiter
,method = method
,indirect.interval = indirect.interval
,indirect.subdivisions = indirect.subdivisions
,indirect.tol = indirect.tol
)$c.i
} else {
c.i <- c.i
}
LCL <- X.bar.bar - c.i * sigma.v
UCL <- X.bar.bar + c.i * sigma.v
if (plot.option == TRUE){
plot(c(1, m), c(LCL, UCL), xaxt = "n", xlab = 'Subgroup', ylab = 'Sample Mean', type = 'n')
axis(side = 1, at = 1:m)
points(1:m, X.bar, type = 'o')
points(c(-1, m + 2), c(LCL, LCL), type = 'l', col = 'red')
points(c(-1, m + 2), c(UCL, UCL), type = 'l', col = 'red')
points(c(-1, m + 2), c(X.bar.bar, X.bar.bar), type = 'l', col = 'blue')
text(round(m * 0.8), UCL, paste('UCL = ', round(UCL, 4)), pos = 1)
text(round(m * 0.8), LCL, paste('LCL = ', round(LCL, 4)), pos = 3)
}
list(CL = X.bar.bar, sigma = sigma.v, PH1.cc = c.i, m = m, nu = nu, LCL = LCL, UCL = UCL, CS = X.bar)
}
####################################################################################################################################################
#Example about table 6.3 (Montgomery, 2013)
####################################################################################################################################################
PH1XBAR.data <- function(){
matrix(c(
74.03 , 74.002 , 74.019 , 73.992 , 74.008 ,
73.995 , 73.992 , 74.001 , 74.011 , 74.004 ,
73.988 , 74.024 , 74.021 , 74.005 , 74.002 ,
74.002 , 73.996 , 73.993 , 74.015 , 74.009 ,
73.992 , 74.007 , 74.015 , 73.989 , 74.014 ,
74.009 , 73.994 , 73.997 , 73.985 , 73.993 ,
73.995 , 74.006 , 73.994 , 74 , 74.005 ,
73.985 , 74.003 , 73.993 , 74.015 , 73.988 ,
74.008 , 73.995 , 74.009 , 74.005 , 74.004 ,
73.998 , 74 , 73.99 , 74.007 , 73.995 ,
73.994 , 73.998 , 73.994 , 73.995 , 73.99 ,
74.004 , 74 , 74.007 , 74 , 73.996 ,
73.983 , 74.002 , 73.998 , 73.997 , 74.012 ,
74.006 , 73.967 , 73.994 , 74 , 73.984 ,
74.012 , 74.014 , 73.998 , 73.999 , 74.007 ,
74 , 73.984 , 74.005 , 73.998 , 73.996 ,
73.994 , 74.012 , 73.986 , 74.005 , 74.007 ,
74.006 , 74.01 , 74.018 , 74.003 , 74 ,
73.984 , 74.002 , 74.003 , 74.005 , 73.997 ,
74 , 74.01 , 74.013 , 74.02 , 74.003 ,
73.982 , 74.001 , 74.015 , 74.005 , 73.996 ,
74.004 , 73.999 , 73.99 , 74.006 , 74.009 ,
74.01 , 73.989 , 73.99 , 74.009 , 74.014 ,
74.015 , 74.008 , 73.993 , 74 , 74.01 ,
73.982 , 73.984 , 73.995 , 74.017 , 74.013
), ncol = 5, byrow = T)
}
####################################################################################################################################################
# Phase II Xbar chart. Please see Yao and Chakraborti (2018)
####################################################################################################################################################
PH2.inner.normal <- function(Z, Y, c.ii, k, nu = k - 1, c4.option = TRUE){
corr.par <- corr.par.f(nu, c4.option)
qn <- Z / sqrt(k) + c.ii / corr.par / sqrt(nu) * sqrt(Y)
pnorm(qn)
}
PH2.CFAR.intgrand <- function(u, v, c.ii, k, nu = k - 1, c4.option = TRUE){
Z <- qnorm(u)
Y <- qchisq(v, nu)
p1 <- PH2.inner.normal(Z, Y, c.ii, k, nu, c4.option)
p2 <- PH2.inner.normal(Z, Y, -c.ii, k, nu, c4.option)
1 - p1 + p2
}
PH2.CARL.intgrand <- function(u, v, c.ii, k, nu = k - 1, c4.option = TRUE) 1 / PH2.CFAR.intgrand(u, v, c.ii, k, nu, c4.option)
####################################################################################################################################################
PH2.get.cc.uc <- function(ARL0, k, nu = k - 1, c4.option = TRUE, interval = c(1, 10), u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
PH2.root.finding.uc <- function(c.ii, k, nu = k - 1, c4.option = TRUE, ARL0, u, v) {
ARL0 - mean(PH2.CARL.intgrand(u, v, c.ii, k, nu, c4.option))
}
k <- k
rt <- uniroot(PH2.root.finding.uc, interval = interval, k = k, nu = nu, c4.option = c4.option, ARL0 = ARL0, u = u, v = v, tol = tol)$root
rt
}
####################################################################################################################################################
PH2.root.finding.EPC <- function(p, k, nu = k - 1, c.ii, ARLb, c4.option = TRUE, u = runif(100000), v = runif(100000)){
1 - p - mean(PH2.CARL.intgrand(u, v, c.ii, k, nu, c4.option) >= ARLb)
}
PH2.get.cc.EPC <- function(p, k, nu = k - 1, eps = 0.1, ARL0 = 370, c4.option = TRUE, interval = c(1, 10), u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
ARLb <- (1 - eps) * ARL0
rt <- uniroot(
PH2.root.finding.EPC,
interval = interval,
p = p,
k = k,
nu = nu,
ARLb = ARLb,
c4.option = c4.option,
u = u,
v = v,
tol = tol
)$root
rt
}
#rnum <- 100000
#u <- runif(rnum)
#v <- runif(rnum)
#debug(PH2.get.c.ii.EPC)
#PH2.get.c.ii.EPC(p = 0.05, k = 100, eps = 0.1, ARL0 = 500, c4.option = TRUE, interval = c(1, 5), u = u, v = v)
PH2.get.ARLb.EPC <- function(p, k, nu = k - 1, c.ii, ARL0 = NULL, c4.option = TRUE, interval = c(1, 1000), u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
rt <- uniroot(
PH2.root.finding.EPC,
interval = interval,
p = p,
k = k,
nu = nu,
c.ii = c.ii,
c4.option = c4.option,
u = u,
v = v,
tol = tol
)$root
if (is.null(ARL0)) {
list(ARLb = rt, eps = NULL, ARL0 = NULL)
} else {
list(ARLb = rt, eps = 1 - rt / ARL0, ARL0 = ARL0)
}
}
#PH2.get.ARLb.EPC(p = 0.1, k = 100, c.ii = 3, c4.option = TRUE, u = u, v = v)
#PH2.get.ARLb.EPC(p = 0.1, k = 100, c.ii = 3, ARL0 = 370, c4.option = TRUE, u = u, v = v)
PH2.get.k.EPC <- function(p, c.ii, eps = 0.1, ARL0 = 370, c4.option = TRUE, interval = c(1000, 3000), model = 'ANOVA-based', n = 10, u = runif(100000), v = runif(100000), tol = .Machine$double.eps^0.25){
if (model == 'ANOVA-based') {
PH2.root.finding.k.EPC <- function(p, k, c.ii, ARLb, c4.option = TRUE, n = 10, u, v){
nu <- k - 1
PH2.root.finding.EPC(p, k, nu, c.ii, ARLb, c4.option, u, v)
}
} else if (model == 'basic') {
PH2.root.finding.k.EPC <- function(p, k, c.ii, ARLb, c4.option = TRUE, n = 10, u, v){
nu <- (n - 1) * k
PH2.root.finding.EPC(p, k, nu, c.ii, ARLb, c4.option, u, v)
}
} else {
stop("need to specify whether it is based on the ANOVA-based model or others")
}
ARLb <- (1 - eps) * ARL0
rt <- uniroot(
PH2.root.finding.k.EPC,
interval = interval,
p = p,
c.ii = c.ii,
ARLb = ARLb,
c4.option = c4.option,
n = n,
u = u,
v = v,
tol = tol
)$root
list(k = rt, PH1.sample.size = rt * n)
}
####################################################################################################################################################
PH2XBAR <- function(
X
,PH1.info = list(X = NULL, mu = NULL, sigma = NULL, k = NULL, n = NULL, model = "ANOVA-based")
,c.ii.info = list(c.ii = NULL, method = 'UC', ARL0 = 370, p = 0.05, eps = 0.1, interval.c.ii.UC = c(1, 3.2), interval.c.ii.EPC = c(1, 10), UC.tol = .Machine$double.eps^0.25, EPC.tol = .Machine$double.eps^0.25)
#,alternative = '2-sided'
,c4.option = TRUE
,plot.option = TRUE
,maxsim = 100000
) {
alternative = '2-sided' #turn off the alternative
k.ii <- dim(X)[1]
X.bar <- rowMeans(X)
if (is.null(PH1.info$X)){
k <- PH1.info$k
if (PH1.info$model == 'ANOVA-based') {
nu <- k - 1
} else if (PH1.info$model == 'basic') {
if (is.null(PH1.info$n)) {
stop('Subgroup size needs to be defined')
} else {
n <- PH1.info$n
nu <- k * (n - 1)
}
} else {
stop('The model in Phase 1 needs to be defined')
}
corr.par <- corr.par.f(nu, c4.option)
X.bar.bar <- PH1.info$mu
sigma.v <- PH1.info$sigma / corr.par
} else {
PH1.chart <- PH1XBAR(X = PH1.info$X, c.i = 1, model = PH1.info$model, plot.option = FALSE)
k <- PH1.chart$k
nu <- PH1.chart$nu
X.bar.bar <- PH1.chart$CL
sigma.v <- PH1.chart$sigma
}
u <- runif(maxsim)
v <- runif(maxsim)
UC <- NULL
EPC <- NULL
n.c.ii <- 0
UC.flag <- 0
EPC.flag <- 0
if (is.null(c.ii.info$c.ii)) {
if (c.ii.info$method == 'both') {
UC.flag <- 1
EPC.flag <- 1
} else if (c.ii.info$method == 'UC') {
UC.flag <- 1
} else if (c.ii.info$method == 'EPC'){
EPC.flag <- 1
} else {
stop("Please check the method of obtaining charting constants in Phase 2")
}
if (UC.flag == 1) {
n.c.ii <- n.c.ii + 1
UC <- PH2.get.cc.uc(
ARL0 = c.ii.info$ARL0
, k = k
, nu = nu
, c4.option = c4.option
, interval = c.ii.info$interval.c.ii.UC
, u = u
, v = v
, tol = c.ii.info$UC.tol
)
}
if (EPC.flag == 1){
n.c.ii <- n.c.ii + 1
EPC <- PH2.get.cc.EPC(
p = c.ii.info$p
, k = k
, nu = nu
, eps = c.ii.info$eps
, ARL0 = c.ii.info$ARL0
, c4.option = c4.option
, interval = c.ii.info$interval.c.ii.EPC
, u = u
, v = v
, tol = c.ii.info$EPC.tol
)
}
c.ii <- list(UC = UC, EPC = EPC)
} else {
n.c.ii <- length(c.ii.info$c.ii)
c.ii <- c.ii.info$c.ii
}
LCL <- X.bar.bar - unlist(c.ii) * sigma.v
UCL <- X.bar.bar + unlist(c.ii) * sigma.v
if (plot.option == TRUE){
plot(c(1, k.ii), c(min(c(LCL, X.bar)), max(c(UCL, X.bar))), xaxt = "n", xlab = 'Subgroup', ylab = 'Sample Mean', type = 'n')
points(1:k.ii, X.bar, type = 'o')
axis(side = 1, at = 1:k.ii)
points(c(-1, k.ii + 2), c(X.bar.bar, X.bar.bar), type = 'l', col = 'blue')
for (ii in 1:n.c.ii){
points(c(-1, k.ii + 2), c(LCL[ii], LCL[ii]), type = 'l', col = 'red')
points(c(-1, k.ii + 2), c(UCL[ii], UCL[ii]), type = 'l', col = 'red')
if (is.null(c.ii.info$c.ii)) {
lab.in <- paste('(', names(c.ii)[ii], ')', sep = '')
} else {
lab.in <- NULL
}
ulab <- paste('UCL', lab.in, ' = ', round(UCL[ii], 4), sep = '')
llab <- paste('LCL', lab.in, ' = ', round(LCL[ii], 4), sep = '')
text(round(k.ii * 0.8), UCL[ii], ulab, pos = 1)
text(round(k.ii * 0.8), LCL[ii], llab, pos = 3)
}
}
list(CL = X.bar.bar, sigma = sigma.v, k = k, nu = nu, PH2.cc = unlist(c.ii), LCL = LCL, UCL = UCL, CS = X.bar)
}
####################################################################################################################################################
PH2XBAR.data <- function(){
matrix(
c(
240, 243, 250, 253, 248,
238, 242, 245, 251, 247,
239, 242, 246, 250, 248,
235, 237, 246, 249, 246,
240, 241, 246, 247, 249,
240, 243, 244, 248, 245,
240, 243, 244, 249, 246,
245, 250, 250, 247, 248,
238, 240, 245, 248, 246,
240, 242, 246, 249, 248,
240, 243, 246, 250, 248,
241, 245, 243, 247, 245,
247, 245, 255, 250, 249,
237, 239, 243, 247, 246,
242, 244, 245, 248, 245,
237, 239, 242, 247, 245,
242, 244, 246, 251, 248,
243, 245, 247, 252, 249,
243, 245, 248, 251, 250,
244, 246, 246, 250, 246,
241, 239, 244, 250, 246,
242, 245, 248, 251, 249,
242, 245, 248, 243, 246,
241, 244, 245, 249, 247,
236, 239, 241, 246, 242,
243, 246, 247, 252, 247,
241, 243, 245, 248, 246,
239, 240, 242, 243, 244,
239, 240, 250, 252, 250,
241, 243, 249, 255, 253
)
,nrow = 30
,ncol = 5
,byrow = TRUE
)
}
#PH2.get.k.EPC(p = 0.05, c.ii = 3, eps = 0.2, ARL0 = 370, c4.option = TRUE, interval = c(100, 400), type = 'ANOVA-based', subgroup.size = 10, u = u, v = v)
#
#m <- dim(x)[1]
#n <- dim(x)[2]
#
#x.bar.bar <- mean(x)
#
#sigma.v <- sqrt(sum(apply(x, 1, var)) / m) / c4.f(m * (n - 1))
#
#k <- get.cc(m, m * (n - 1), 0.05)$K
#
#LCL <- x.bar.bar - k * sigma.v / sqrt(n)
#UCL <- x.bar.bar + k * sigma.v / sqrt(n)
#
#
#x.bar <- rowMeans(x)
#
#plot(c(1, m), c(LCL * 49997/50000, UCL * 50003/50000), xaxt = "n", xlab = 'Subgroup', ylab = 'Sample Means', type = 'n')
#
#axis(side = 1, at = 1:25)
#
#points(1:m, x.bar, type = 'o')
#points(c(-1, m + 2), c(LCL, LCL), type = 'l', col = 'red')
#points(c(-1, m + 2), c(UCL, UCL), type = 'l', col = 'red')
#points(c(-1, m + 2), c(x.bar.bar, x.bar.bar), type = 'l', col = 'blue')
#text(20, UCL * 50001/50000, paste('UCL = ', round(UCL, 4)))
#text(20, LCL * 49999/50000, paste('LCL = ', round(LCL, 4)))
#
#pos <- c(rep(c(3, 1), 5), 1, 3, 3, 1, 3, 1, 1, 3, 1, 3, 1, 3, 1, 3, 1) #rep(c(1, 3), 8))[-26]
#
#text(1:m, x.bar, round(x.bar, 4), cex = 0.7, pos = pos)
#
#
####################################################################################################################################################
#Example about how to get K and L by FAP0
####################################################################################################################################################
#get.cc(
# 20
# ,80
# ,FAP = 0.1
# ,off.diag = NULL
# ,alternative = '2-sided'
# ,maxiter = 10000
# ,method = 'direct'
#)
#get.cc(
# 20
# ,80
# ,FAP = 0.1
# ,off.diag = NULL
# ,alternative = '2-sided'
# ,maxiter = 10000
# ,method = 'indirect'
#)
####################################################################################################################################################
#Example about how to get FAP by K
####################################################################################################################################################
#
#k <- c(
#22.84794
#,10.52163
#,7.386245
#,6.057107
#,5.340985
#,4.898076
#,4.598652
#,4.383316
#,3.842428
#,3.619871
#,3.498807
#,3.422752
#,3.280788
#,3.182391
#,3.151007
#,3.135567
#,3.126384
#,3.120294
#
#
#
#
#
#)
#
#m.seq <- c(
# 3
# ,4
# ,5
# ,6
# ,7
# ,8
# ,9
# ,10
# ,15
# ,20
# ,25
# ,30
# ,50
# ,100
# ,150
# ,200
# ,250
# ,300
#
#)
#
#
#i <- 0
#
#record <- rep(NA, 18)
#
#for (m in m.seq){
#
# i <- i + 1
#
# nu <- m - 1
#
# record[i] <- get.FAP0(k[i], m, nu)
#
#}
#
#as.matrix(record)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.