R/PCC_function_Gamma.R
In BayesianSPCM/BSPCM: Bayesian Statistical Process Control and Monitoring
Defines functions
PCC_Gamma
Documented in
PCC_Gamma
# Function for running PCC Gamma with unknown rate parameter
PCC_Gamma <- function(data = NULL, historical_data = NULL, al = NULL,
c0 = 0, d0 = 0, alpha_0 = NULL,
ARL_0 = 370.4, FAP = NULL, FIR = FALSE, fFIR = .99,
aFIR = 1/8, summary_list = TRUE, PCC_PLOT = TRUE, PriorPosterior_PLOT = FALSE,
historical_data_PLOT = FALSE, pdf_report = FALSE,
xlab = "Observation Order", ylab = "Quality characteristic Values",
main="PCC Gamma with rate unknown"){
### Initial checks before procceeding to the main body of function
### Mainly this chunk of code will correspond to invalid general input before running stuff
# 'data' (i) not defined (ii) not in vector (iii) contain non-numeric value
if ( is.null(data) ) {
stop("'data' have not been defined")
} else { if ( any(!is.numeric((unlist(data)))) ) stop("Invalid 'data' input")
if ( !is.vector(data) ) stop("'data' must be in vector form")
}
# 'historical_data' (i) not in vector (ii) contain non-numeric value
if ( !is.null(historical_data) ) {
if ( any(!is.numeric((unlist(historical_data)))) ) stop("Invalid 'historical_data' input")
if ( !is.vector(data) ) stop("'historical data' must be in vector form")
}
# 'ARL_0' (i) non-numeric (ii) negative
if( !missing(ARL_0) ) {
if ( length(unlist(ARL_0))>1 ) { message("More than one value for 'ARL_0', the first one will only be used")
if ( !is.numeric(ARL_0[1]) | ARL_0<=0 ) { stop("Invalid 'ARL_0' value") } else { ARL_0 <- ARL_0[1] }
} else { if ( !is.numeric(ARL_0) | ARL_0<=0 ) { stop("Invalid 'ARL_0' value") } }
}
# 'FAP' (i) non-numeric (ii) negative
if (!missing(FAP)){
if ( length(unlist(FAP))>1 ) { message("More than one value for 'FAP', the first one will only be used")
if ( !is.numeric(FAP[1]) | FAP<=0 | FAP>=1 ) { stop("Invalid 'FAP' value") } else { FAP <- FAP[1] }
} else { if ( !is.numeric(FAP) | FAP<=0 | FAP>=1 ) { stop("Invalid 'FAP' value") } }
}
# 'FIR' (i) logical (ii) fFIR - aFIR conditions
if ( length(unlist(FIR))>1 ) {
message("More than one value for 'FIR', the first one will only be used")
if ( !is.logical(FIR[1]) ) { stop("Invalid 'FIR' value ; 'FIR' must be logical") } else { FIR <- FIR[1] }
} else {
if ( !is.logical(FIR) ) { stop("Invalid 'FIR' value ; 'FIR' must be logical") }
}
# fFIR - aFIR conditions if FIR
if ( FIR ) {
if ( !missing(fFIR) ) {
if ( length(unlist(fFIR))>1 ) {
message("More than one value for 'fFIR', the first one will only be used")
if ( !is.numeric(fFIR[1]) | fFIR[1]<=0 | fFIR[1]>=1 ) {
stop("Invalid 'fFIR' value")
} else { fFIR <- fFIR[1] }
} else {
if ( !is.numeric(fFIR) | fFIR<=0 | fFIR>=1 ) {
stop("Invalid 'fFIR' value")
}
}
}
if ( !missing(aFIR) ) {
if ( length(unlist(aFIR))>1 ) {
message("More than one value for 'aFIR', the first one will only be used")
if ( !is.numeric(aFIR[1]) | aFIR[1]<=0 ) {
stop("Invalid 'aFIR' value")
} else { aFIR <- aFIR[1] }
} else {
if ( !is.numeric(aFIR) | aFIR<=0 ) {
stop("Invalid 'aFIR' value")
}
}
}
}
### Setting the False Alarm Probability & False Alarm Rate based on the Sidak correction
# data length
N <- length(data)
# If both ARL_0 and FAP chosen
if ( !is.null(ARL_0) & !is.null(FAP) ) {
message("Both ARL_0 and FAP are defined as input, so ARL_0 is used by default. \nIn order to use FAP instead, set ARL_0 = NULL")
FAR <- 1/ARL_0
# If only FAP is chosen
} else if ( is.null(ARL_0) & !is.null(FAP) ) {
FAR <- 1-(1-FAP)^(1/(N-1))
# If only ARL0 is chosen
} else if ( !is.null(ARL_0) & is.null(FAP) ){
FAR <- 1/ARL_0
}
# If FIR PCC is chosen - default value for f=0.99
if ( FIR ) {
tf <- 1:N
Afir <- c(( 1- (1-fFIR)^(1+aFIR*(tf-1)) ) )
FAR <- 1-(1-FAR)*Afir
}
###############################################################
###############################################################
## START (1) Only this bit changes from function to function ##
## In some cases 'data' and 'historical_data' restrictions ##
## change as well at the beginning of the function ##
###############################################################
###############################################################
# Likelihood al input (i) more than one value for parameters (ii) non-numeric input
if ( is.null(al) ) {
stop("Likelihood shape parameter 'al' has not been defined")
} else {
if ( length(unlist(al))>1 ) { message("More than one value for 'al', the first one will only be used")
if ( !is.numeric(al) | al<0 ) { stop("Invalid 'al' value") } else { al <- al[1] }
} else { if ( !is.numeric(al) | al<0 ) { stop("Invalid 'al' value") } }
}
# Prior parameter input (i) more than one value for parameters (ii) non-numeric input
if( !missing(c0) ) {
if ( length(unlist(c0))>1 ) { message("More than one value for 'c0', the first one will only be used")
if ( !is.numeric(c0) | c0<0 ) { stop("Invalid 'c0' value") } else { c0 <- c0[1] }
} else { if ( !is.numeric(c0) | c0<0 ) { stop("Invalid 'c0' value") } }
}
if( !missing(d0) ) {
if ( length(unlist(d0))>1 ) { message("More than one value for 'd0', the first one will only be used")
if ( !is.numeric(d0) | d0<0 ) { stop("Invalid 'd0' value") } else { d0 <- d0[1] }
} else { if ( !is.numeric(d0) | d0<0 ) { stop("Invalid 'd0' value") } }
}
### Main body of function - PCC illustration - USING FAR (or FAP equivelantly)
## Histotic data and processing
if ( !is.null(historical_data) ){
N_historicaldata <- length(historical_data)
# Check about alpha_0
# If no chosen value for alpha_0 use default setting
if (is.null(alpha_0)) { alpha_0 <-1/N_historicaldata
} else {
if ( length(unlist(alpha_0))>1 ) {
message("More than one value for 'alpha_0', the first one will only be used")
if ( !is.numeric(alpha_0) | alpha_0<0 | alpha_0>1) { stop("Invalid 'alpha_0' value")
} else { if ( !is.numeric(alpha_0) | alpha_0<0 | alpha_0>1 ) { stop("Invalid 'alpha_0' value") } }
}
}
# Process historical data
# Power Prior parameters
c0_PowerP <- c0 + a*alpha_0*N_historicaldata
d0_PowerP <- d0 + alpha_0*sum(historical_data)
# Keep similar notation as input
c0 <- c0_PowerP ; d0 <- d0_PowerP
}
### PCC implementation
# Sum of observations
dataSum <- cumsum(data)[seq(1, length(data))]
# Sum of squared observations
dataSumSquared <- cumsum(data^2)[seq(1, length(data))]
# Posterior parameters
c0_Post <- c0 + al*(1:N)
d0_Post <- d0 + dataSum
# HPRD - Gamma likelihood ; rate unknown
HPRD_likGamRtU <- function(far, A, Cp, Dp){
f <- function(x) {extraDistr::qbetapr( 1-far+x, shape1 = A, shape2 = Cp, scale = Dp ) - extraDistr::qbetapr( x, shape1 = A, shape2 = Cp, scale = Dp )}
out <- optimize( f, c(0,far), tol = .Machine$double.eps )
c( extraDistr::qbetapr(out$minimum, shape1 = A, shape2 = Cp, scale = Dp), extraDistr::qbetapr(1-far+out$minimum, shape1 = A, shape2 = Cp, scale = Dp) )
}
# Control limits
if (!FIR) { CL <- t( mapply( function(A0, C0, D0, FD=FAR) { HPRD_likGamRtU(far=FD, A=A0, Cp=C0, Dp=D0) }
, A0=al, C0=c0_Post, D0=d0_Post) )
} else { CL <- t( mapply( function(A0, C0, D0, FD) { HPRD_likGamRtU(far=FD, A=A0, Cp=C0, Dp=D0) }
, A0=al, C0=c0_Post, D0=d0_Post, FD=FAR) )
}
CL <- rbind( c(NA, NA), CL )
CL <- CL[-nrow(CL), ]
if ( PriorPosterior_PLOT ) {
PrPostPLOT <-
ggplot2::ggplot(data.frame(x = c(
ifelse( c0 == 0 & d0 == 0, .15,
min( qgamma(1 - .999, shape = c0, rate = d0), qgamma(1 - .999, shape = c0_Post[N], rate = d0_Post[N]) ) ),
ifelse( c0 == 0 & d0 == 0, qgamma(.999, shape = c0_Post[N], rate = d0_Post[N]) + c0_Post[N]/d0_Post[N],
if (c0 > d0) { min( qgamma(.999, shape = c0, rate = d0), qgamma(.999, shape = c0_Post[N], rate = d0_Post[N]) ) + c0/d0
} else { max( qgamma(.999, shape = c0, rate = d0), qgamma(.999, shape = c0_Post[N], rate = d0_Post[N]) ) + c0/d0 } )
)
),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = dgamma, args = list(shape = c0_Post[N], rate = d0_Post[N]),
aes(colour = "Posterior", linetype = "Posterior"), size = 1) +
{if(c0 == 0 & d0 == 0) { ggplot2::stat_function(fun = function(x) { 1/(x + abs(log(c0_Post[N]/d0_Post[N]))) }, ggplot2::aes(colour = "Prior", linetype = "Prior"), size = 1)
} else { ggplot2::stat_function(fun = dgamma, args = list(shape = c0, rate = d0), ggplot2::aes(colour = "Prior", linetype = "Prior"), size = 1) } } +
ggplot2::scale_x_continuous(name = "") +
ggplot2::scale_y_continuous(name = "Density") +
ggplot2::scale_linetype_manual(values = c("solid", "dashed"), guide = FALSE) +
ggplot2::scale_colour_manual(values = c("#3CB371", "#FF4500"),
labels = c( bquote("Prior: Gamma(" ~ theta ~ "|" ~ .(round(c0, digits = 1)) ~ ", " ~ .(round(d0, digits = 1)) ~ ")" ),
bquote("Posterior: Gamma(" ~ theta ~ "|" ~ .(round(c0_Post[N], digits = 1)) ~ ", " ~ .(round(d0_Post[N], digits = 1)) ~ ")" )),
guide = guide_legend(override.aes = list( color = c("#FF4500", "#3CB371"),
linetype = c("dashed", "solid"),
size = c(.5, .5)), title = NULL)) +
ggplot2::ggtitle(expression(atop("PCC Gamma likelihood - unknown rate"~theta, "Prior/Posterior distribution"))) +
{if(c0 == 0 & d0 == 0) {
ggplot2::geom_point( ggplot2::aes(x = c0_Post[N]/d0_Post[N], y = 0), color = "#3CB371", show.legend = FALSE, shape = 4, size = 3, stroke = 1.5, na.rm = TRUE )
} else { ggplot2::geom_point( ggplot2::aes(x = c(c0/d0, c0_Post[N]/d0_Post[N]), y = c(0, 0)), color = c("#FF4500", "#3CB371"), show.legend = FALSE, shape = 4, size = 3, stroke = 1.5, na.rm = TRUE ) } } +
{if(c0 == 0 & d0 == 0) {
ggplot2::annotate("text", x = c0_Post[N]/d0_Post[N], y = 0, label = paste(expression(mu[post])), color = "#3CB371", size = 6, parse = TRUE, vjust = 1.25)
} else { ggplot2::annotate("text", x = c(c0/d0, c0_Post[N]/d0_Post[N]), y = c(0, 0), label = paste(expression(mu[prior], mu[post])),
color = c("#FF4500", "#3CB371"), size = 6, parse = TRUE, vjust = 1.25) } } +
ggplot2::theme(legend.position = "bottom",
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 18, hjust = 0.5),
text = element_text(size = 15),
axis.text.x = element_text(colour="black", size = 12),
axis.text.y = element_text(colour="black", size = 12))
print(PrPostPLOT)
}
####################################################################
####################################################################
## END (1) Only the above bit changes from function to function ##
####################################################################
####################################################################
## Output
{ # Construction of 'In' and 'Out' of control column for return results
States <- rep("", times=N)
States[ifelse(data < CL[, 1], TRUE, FALSE)] <- "Alarm (LL)" ; States[ifelse(data > CL[, 2], TRUE, FALSE)] <- "Alarm (UL)"
# Return results
PCC_summary <- data.frame( data=data, HPrD_LL=CL[, 1], HPrD_UL=CL[, 2], Alarms=States ) }
## Dynamic recalculation of PCC plot's y axis
# PCC y axis limits allowance
Ratio <- (PCC_summary$HPrD_UL-PCC_summary$HPrD_LL)/min(PCC_summary$HPrD_UL-PCC_summary$HPrD_LL, na.rm = T)
# Y axis limits
AdjustedYlim <- c(min(PCC_summary$data, PCC_summary$HPrD_LL[which(Ratio<=2.5)], na.rm=T),
max(PCC_summary$data, PCC_summary$HPrD_UL[which(Ratio<=2.5)], na.rm=T))
### Output of function
## PCC plot
if ( PCC_PLOT ) {
# Creation of PCC plot
PCC_PlotSummary <- cbind(Indices=1:N, PCC_summary)
PCC <- ggplot2::ggplot(PCC_PlotSummary, ggplot2::aes(Indices, data)) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=data), na.rm = TRUE) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_UL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_LL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_ribbon(ggplot2::aes(x=Indices, ymin=HPrD_UL, ymax=HPrD_LL, fill=TRUE), alpha=0.25, show.legend=FALSE) +
ggplot2::scale_fill_manual(values=c("TRUE"="green")) +
ggplot2::geom_point(ggplot2::aes(group=Indices, color=as.factor(Alarms), stroke = 1.5), show.legend=FALSE, na.rm = TRUE) +
ggplot2::scale_color_manual(values=c("black", "red", "red"), na.value = "black") +
ggplot2::coord_cartesian(ylim = AdjustedYlim) +
ggplot2::labs(title = main, x = xlab, y = ylab) +
ggplot2::theme(legend.position = "top",
legend.title = ggplot2::element_blank(),
axis.line = ggplot2::element_line(colour = "black", size = 0.5, linetype = "solid"),
panel.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5)
)
# Creation of PCC plot if historical data are chosen to be on the plot
if ( !is.null(historical_data) & historical_data_PLOT ) {
PCC_summary_historicaldata <- data.frame( data=c(historical_data, data), HPrD_LL=c(rep(NA, times=N_historicaldata), CL[, 1]),
HPrD_UL=c(rep(NA, times=N_historicaldata), CL[, 2]), Alarms=c(rep("", times=N_historicaldata), States) )
PCC_PlotSummary <- cbind(Indices=c(-N_historicaldata:(-1), 1:N), TypeOfdata=c(rep("Historical", times=N_historicaldata), rep("Current", times=N)), PCC_summary_historicaldata )
PCC_historical <- ggplot2::ggplot(PCC_PlotSummary, ggplot2::aes(Indices, data)) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=data, linetype = as.factor(TypeOfdata)), na.rm = TRUE) +
ggplot2::geom_segment(ggplot2::aes(x = 0, y = min(HPrD_LL, na.rm=TRUE), xend = 0, yend = max(HPrD_UL, na.rm=TRUE))) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_UL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_LL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_ribbon(ggplot2::aes(x=Indices, ymin=HPrD_UL, ymax=HPrD_LL, fill=TRUE), alpha=0.25, show.legend=FALSE) +
ggplot2::scale_fill_manual(values=c("TRUE"="green")) +
ggplot2::geom_point(ggplot2::aes(group=Indices, shape=as.factor(TypeOfdata), color=as.factor(Alarms), stroke = 1.5), show.legend=FALSE, na.rm = TRUE) +
ggplot2::scale_color_manual(values=c("black", "red", "red"), na.value = "black") +
ggplot2::scale_linetype_manual(values=c("Historical"="dotted", "Current"="solid")) +
ggplot2::scale_shape_manual(values=c("Historical"=1, "Current"=19)) +
ggplot2::coord_cartesian(ylim = AdjustedYlim) +
ggplot2::labs(title = main, x = xlab, y = ylab) +
ggplot2::theme(legend.position = "top",
legend.title = ggplot2::element_blank(),
axis.line = ggplot2::element_line(colour = "black", size = 0.5, linetype = "solid"),
panel.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5)
)
print(PCC_historical)
} else { print(PCC) }
}
# List of results
if ( summary_list ) { print(PCC_summary) }
# List of results return in pdf
if ( pdf_report ) {
# save pdf
Cairo::CairoPDF(
paste0( "PCC_results_", paste0( unlist(strsplit(date(), " "))[c(1,2,3,5)], collapse = "_" ), "_",
paste0( unlist(strsplit( unlist(strsplit(date(), " "))[4], ":" )), collapse = "." ),
".pdf" ),
height = 8.264, width = 11.694)
# PCC plot on pdf
print(PCC)
# Prior Posterior plot on pdf
print(PrPostPLOT)
# Results matrix on pdf
# Chunk of code to split results matrix to different pages - Set a default number based on pdf height/width
NRowsPerPage <- 25
if(NRowsPerPage > nrow(PCC_summary)){ FloatingRow <- nrow(PCC_summary) } else { FloatingRow <- NRowsPerPage }
sapply(1:ceiling(nrow(PCC_summary)/NRowsPerPage), function(index) {
if (index==1) { StartingRow <<- 1 }
grid::grid.newpage()
gridExtra::grid.table(PCC_summary[StartingRow:FloatingRow, ])
StartingRow <<- FloatingRow + 1
if( sum(NRowsPerPage, FloatingRow) < nrow(PCC_summary)){ FloatingRow <<- NRowsPerPage + FloatingRow } else { FloatingRow <<- nrow(PCC_summary) }
})
grDevices::dev.off()
}
}
#out <- rgamma(n = 30, shape = 4, rate = 4)
#out[15] <- out[15] + 1.5
#PCC_Gamma(out, al = 4, PriorPosterior_PLOT = TRUE, PCC_PLOT = FALSE)
BayesianSPCM/BSPCM documentation built on Dec. 31, 2020, 9:53 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 PCC_Gamma
Documented in PCC_Gamma
# Function for running PCC Gamma with unknown rate parameter
PCC_Gamma <- function(data = NULL, historical_data = NULL, al = NULL,
c0 = 0, d0 = 0, alpha_0 = NULL,
ARL_0 = 370.4, FAP = NULL, FIR = FALSE, fFIR = .99,
aFIR = 1/8, summary_list = TRUE, PCC_PLOT = TRUE, PriorPosterior_PLOT = FALSE,
historical_data_PLOT = FALSE, pdf_report = FALSE,
xlab = "Observation Order", ylab = "Quality characteristic Values",
main="PCC Gamma with rate unknown"){
### Initial checks before procceeding to the main body of function
### Mainly this chunk of code will correspond to invalid general input before running stuff
# 'data' (i) not defined (ii) not in vector (iii) contain non-numeric value
if ( is.null(data) ) {
stop("'data' have not been defined")
} else { if ( any(!is.numeric((unlist(data)))) ) stop("Invalid 'data' input")
if ( !is.vector(data) ) stop("'data' must be in vector form")
}
# 'historical_data' (i) not in vector (ii) contain non-numeric value
if ( !is.null(historical_data) ) {
if ( any(!is.numeric((unlist(historical_data)))) ) stop("Invalid 'historical_data' input")
if ( !is.vector(data) ) stop("'historical data' must be in vector form")
}
# 'ARL_0' (i) non-numeric (ii) negative
if( !missing(ARL_0) ) {
if ( length(unlist(ARL_0))>1 ) { message("More than one value for 'ARL_0', the first one will only be used")
if ( !is.numeric(ARL_0[1]) | ARL_0<=0 ) { stop("Invalid 'ARL_0' value") } else { ARL_0 <- ARL_0[1] }
} else { if ( !is.numeric(ARL_0) | ARL_0<=0 ) { stop("Invalid 'ARL_0' value") } }
}
# 'FAP' (i) non-numeric (ii) negative
if (!missing(FAP)){
if ( length(unlist(FAP))>1 ) { message("More than one value for 'FAP', the first one will only be used")
if ( !is.numeric(FAP[1]) | FAP<=0 | FAP>=1 ) { stop("Invalid 'FAP' value") } else { FAP <- FAP[1] }
} else { if ( !is.numeric(FAP) | FAP<=0 | FAP>=1 ) { stop("Invalid 'FAP' value") } }
}
# 'FIR' (i) logical (ii) fFIR - aFIR conditions
if ( length(unlist(FIR))>1 ) {
message("More than one value for 'FIR', the first one will only be used")
if ( !is.logical(FIR[1]) ) { stop("Invalid 'FIR' value ; 'FIR' must be logical") } else { FIR <- FIR[1] }
} else {
if ( !is.logical(FIR) ) { stop("Invalid 'FIR' value ; 'FIR' must be logical") }
}
# fFIR - aFIR conditions if FIR
if ( FIR ) {
if ( !missing(fFIR) ) {
if ( length(unlist(fFIR))>1 ) {
message("More than one value for 'fFIR', the first one will only be used")
if ( !is.numeric(fFIR[1]) | fFIR[1]<=0 | fFIR[1]>=1 ) {
stop("Invalid 'fFIR' value")
} else { fFIR <- fFIR[1] }
} else {
if ( !is.numeric(fFIR) | fFIR<=0 | fFIR>=1 ) {
stop("Invalid 'fFIR' value")
}
}
}
if ( !missing(aFIR) ) {
if ( length(unlist(aFIR))>1 ) {
message("More than one value for 'aFIR', the first one will only be used")
if ( !is.numeric(aFIR[1]) | aFIR[1]<=0 ) {
stop("Invalid 'aFIR' value")
} else { aFIR <- aFIR[1] }
} else {
if ( !is.numeric(aFIR) | aFIR<=0 ) {
stop("Invalid 'aFIR' value")
}
}
}
}
### Setting the False Alarm Probability & False Alarm Rate based on the Sidak correction
# data length
N <- length(data)
# If both ARL_0 and FAP chosen
if ( !is.null(ARL_0) & !is.null(FAP) ) {
message("Both ARL_0 and FAP are defined as input, so ARL_0 is used by default. \nIn order to use FAP instead, set ARL_0 = NULL")
FAR <- 1/ARL_0
# If only FAP is chosen
} else if ( is.null(ARL_0) & !is.null(FAP) ) {
FAR <- 1-(1-FAP)^(1/(N-1))
# If only ARL0 is chosen
} else if ( !is.null(ARL_0) & is.null(FAP) ){
FAR <- 1/ARL_0
}
# If FIR PCC is chosen - default value for f=0.99
if ( FIR ) {
tf <- 1:N
Afir <- c(( 1- (1-fFIR)^(1+aFIR*(tf-1)) ) )
FAR <- 1-(1-FAR)*Afir
}
###############################################################
###############################################################
## START (1) Only this bit changes from function to function ##
## In some cases 'data' and 'historical_data' restrictions ##
## change as well at the beginning of the function ##
###############################################################
###############################################################
# Likelihood al input (i) more than one value for parameters (ii) non-numeric input
if ( is.null(al) ) {
stop("Likelihood shape parameter 'al' has not been defined")
} else {
if ( length(unlist(al))>1 ) { message("More than one value for 'al', the first one will only be used")
if ( !is.numeric(al) | al<0 ) { stop("Invalid 'al' value") } else { al <- al[1] }
} else { if ( !is.numeric(al) | al<0 ) { stop("Invalid 'al' value") } }
}
# Prior parameter input (i) more than one value for parameters (ii) non-numeric input
if( !missing(c0) ) {
if ( length(unlist(c0))>1 ) { message("More than one value for 'c0', the first one will only be used")
if ( !is.numeric(c0) | c0<0 ) { stop("Invalid 'c0' value") } else { c0 <- c0[1] }
} else { if ( !is.numeric(c0) | c0<0 ) { stop("Invalid 'c0' value") } }
}
if( !missing(d0) ) {
if ( length(unlist(d0))>1 ) { message("More than one value for 'd0', the first one will only be used")
if ( !is.numeric(d0) | d0<0 ) { stop("Invalid 'd0' value") } else { d0 <- d0[1] }
} else { if ( !is.numeric(d0) | d0<0 ) { stop("Invalid 'd0' value") } }
}
### Main body of function - PCC illustration - USING FAR (or FAP equivelantly)
## Histotic data and processing
if ( !is.null(historical_data) ){
N_historicaldata <- length(historical_data)
# Check about alpha_0
# If no chosen value for alpha_0 use default setting
if (is.null(alpha_0)) { alpha_0 <-1/N_historicaldata
} else {
if ( length(unlist(alpha_0))>1 ) {
message("More than one value for 'alpha_0', the first one will only be used")
if ( !is.numeric(alpha_0) | alpha_0<0 | alpha_0>1) { stop("Invalid 'alpha_0' value")
} else { if ( !is.numeric(alpha_0) | alpha_0<0 | alpha_0>1 ) { stop("Invalid 'alpha_0' value") } }
}
}
# Process historical data
# Power Prior parameters
c0_PowerP <- c0 + a*alpha_0*N_historicaldata
d0_PowerP <- d0 + alpha_0*sum(historical_data)
# Keep similar notation as input
c0 <- c0_PowerP ; d0 <- d0_PowerP
}
### PCC implementation
# Sum of observations
dataSum <- cumsum(data)[seq(1, length(data))]
# Sum of squared observations
dataSumSquared <- cumsum(data^2)[seq(1, length(data))]
# Posterior parameters
c0_Post <- c0 + al*(1:N)
d0_Post <- d0 + dataSum
# HPRD - Gamma likelihood ; rate unknown
HPRD_likGamRtU <- function(far, A, Cp, Dp){
f <- function(x) {extraDistr::qbetapr( 1-far+x, shape1 = A, shape2 = Cp, scale = Dp ) - extraDistr::qbetapr( x, shape1 = A, shape2 = Cp, scale = Dp )}
out <- optimize( f, c(0,far), tol = .Machine$double.eps )
c( extraDistr::qbetapr(out$minimum, shape1 = A, shape2 = Cp, scale = Dp), extraDistr::qbetapr(1-far+out$minimum, shape1 = A, shape2 = Cp, scale = Dp) )
}
# Control limits
if (!FIR) { CL <- t( mapply( function(A0, C0, D0, FD=FAR) { HPRD_likGamRtU(far=FD, A=A0, Cp=C0, Dp=D0) }
, A0=al, C0=c0_Post, D0=d0_Post) )
} else { CL <- t( mapply( function(A0, C0, D0, FD) { HPRD_likGamRtU(far=FD, A=A0, Cp=C0, Dp=D0) }
, A0=al, C0=c0_Post, D0=d0_Post, FD=FAR) )
}
CL <- rbind( c(NA, NA), CL )
CL <- CL[-nrow(CL), ]
if ( PriorPosterior_PLOT ) {
PrPostPLOT <-
ggplot2::ggplot(data.frame(x = c(
ifelse( c0 == 0 & d0 == 0, .15,
min( qgamma(1 - .999, shape = c0, rate = d0), qgamma(1 - .999, shape = c0_Post[N], rate = d0_Post[N]) ) ),
ifelse( c0 == 0 & d0 == 0, qgamma(.999, shape = c0_Post[N], rate = d0_Post[N]) + c0_Post[N]/d0_Post[N],
if (c0 > d0) { min( qgamma(.999, shape = c0, rate = d0), qgamma(.999, shape = c0_Post[N], rate = d0_Post[N]) ) + c0/d0
} else { max( qgamma(.999, shape = c0, rate = d0), qgamma(.999, shape = c0_Post[N], rate = d0_Post[N]) ) + c0/d0 } )
)
),
ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = dgamma, args = list(shape = c0_Post[N], rate = d0_Post[N]),
aes(colour = "Posterior", linetype = "Posterior"), size = 1) +
{if(c0 == 0 & d0 == 0) { ggplot2::stat_function(fun = function(x) { 1/(x + abs(log(c0_Post[N]/d0_Post[N]))) }, ggplot2::aes(colour = "Prior", linetype = "Prior"), size = 1)
} else { ggplot2::stat_function(fun = dgamma, args = list(shape = c0, rate = d0), ggplot2::aes(colour = "Prior", linetype = "Prior"), size = 1) } } +
ggplot2::scale_x_continuous(name = "") +
ggplot2::scale_y_continuous(name = "Density") +
ggplot2::scale_linetype_manual(values = c("solid", "dashed"), guide = FALSE) +
ggplot2::scale_colour_manual(values = c("#3CB371", "#FF4500"),
labels = c( bquote("Prior: Gamma(" ~ theta ~ "|" ~ .(round(c0, digits = 1)) ~ ", " ~ .(round(d0, digits = 1)) ~ ")" ),
bquote("Posterior: Gamma(" ~ theta ~ "|" ~ .(round(c0_Post[N], digits = 1)) ~ ", " ~ .(round(d0_Post[N], digits = 1)) ~ ")" )),
guide = guide_legend(override.aes = list( color = c("#FF4500", "#3CB371"),
linetype = c("dashed", "solid"),
size = c(.5, .5)), title = NULL)) +
ggplot2::ggtitle(expression(atop("PCC Gamma likelihood - unknown rate"~theta, "Prior/Posterior distribution"))) +
{if(c0 == 0 & d0 == 0) {
ggplot2::geom_point( ggplot2::aes(x = c0_Post[N]/d0_Post[N], y = 0), color = "#3CB371", show.legend = FALSE, shape = 4, size = 3, stroke = 1.5, na.rm = TRUE )
} else { ggplot2::geom_point( ggplot2::aes(x = c(c0/d0, c0_Post[N]/d0_Post[N]), y = c(0, 0)), color = c("#FF4500", "#3CB371"), show.legend = FALSE, shape = 4, size = 3, stroke = 1.5, na.rm = TRUE ) } } +
{if(c0 == 0 & d0 == 0) {
ggplot2::annotate("text", x = c0_Post[N]/d0_Post[N], y = 0, label = paste(expression(mu[post])), color = "#3CB371", size = 6, parse = TRUE, vjust = 1.25)
} else { ggplot2::annotate("text", x = c(c0/d0, c0_Post[N]/d0_Post[N]), y = c(0, 0), label = paste(expression(mu[prior], mu[post])),
color = c("#FF4500", "#3CB371"), size = 6, parse = TRUE, vjust = 1.25) } } +
ggplot2::theme(legend.position = "bottom",
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 18, hjust = 0.5),
text = element_text(size = 15),
axis.text.x = element_text(colour="black", size = 12),
axis.text.y = element_text(colour="black", size = 12))
print(PrPostPLOT)
}
####################################################################
####################################################################
## END (1) Only the above bit changes from function to function ##
####################################################################
####################################################################
## Output
{ # Construction of 'In' and 'Out' of control column for return results
States <- rep("", times=N)
States[ifelse(data < CL[, 1], TRUE, FALSE)] <- "Alarm (LL)" ; States[ifelse(data > CL[, 2], TRUE, FALSE)] <- "Alarm (UL)"
# Return results
PCC_summary <- data.frame( data=data, HPrD_LL=CL[, 1], HPrD_UL=CL[, 2], Alarms=States ) }
## Dynamic recalculation of PCC plot's y axis
# PCC y axis limits allowance
Ratio <- (PCC_summary$HPrD_UL-PCC_summary$HPrD_LL)/min(PCC_summary$HPrD_UL-PCC_summary$HPrD_LL, na.rm = T)
# Y axis limits
AdjustedYlim <- c(min(PCC_summary$data, PCC_summary$HPrD_LL[which(Ratio<=2.5)], na.rm=T),
max(PCC_summary$data, PCC_summary$HPrD_UL[which(Ratio<=2.5)], na.rm=T))
### Output of function
## PCC plot
if ( PCC_PLOT ) {
# Creation of PCC plot
PCC_PlotSummary <- cbind(Indices=1:N, PCC_summary)
PCC <- ggplot2::ggplot(PCC_PlotSummary, ggplot2::aes(Indices, data)) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=data), na.rm = TRUE) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_UL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_LL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_ribbon(ggplot2::aes(x=Indices, ymin=HPrD_UL, ymax=HPrD_LL, fill=TRUE), alpha=0.25, show.legend=FALSE) +
ggplot2::scale_fill_manual(values=c("TRUE"="green")) +
ggplot2::geom_point(ggplot2::aes(group=Indices, color=as.factor(Alarms), stroke = 1.5), show.legend=FALSE, na.rm = TRUE) +
ggplot2::scale_color_manual(values=c("black", "red", "red"), na.value = "black") +
ggplot2::coord_cartesian(ylim = AdjustedYlim) +
ggplot2::labs(title = main, x = xlab, y = ylab) +
ggplot2::theme(legend.position = "top",
legend.title = ggplot2::element_blank(),
axis.line = ggplot2::element_line(colour = "black", size = 0.5, linetype = "solid"),
panel.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5)
)
# Creation of PCC plot if historical data are chosen to be on the plot
if ( !is.null(historical_data) & historical_data_PLOT ) {
PCC_summary_historicaldata <- data.frame( data=c(historical_data, data), HPrD_LL=c(rep(NA, times=N_historicaldata), CL[, 1]),
HPrD_UL=c(rep(NA, times=N_historicaldata), CL[, 2]), Alarms=c(rep("", times=N_historicaldata), States) )
PCC_PlotSummary <- cbind(Indices=c(-N_historicaldata:(-1), 1:N), TypeOfdata=c(rep("Historical", times=N_historicaldata), rep("Current", times=N)), PCC_summary_historicaldata )
PCC_historical <- ggplot2::ggplot(PCC_PlotSummary, ggplot2::aes(Indices, data)) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=data, linetype = as.factor(TypeOfdata)), na.rm = TRUE) +
ggplot2::geom_segment(ggplot2::aes(x = 0, y = min(HPrD_LL, na.rm=TRUE), xend = 0, yend = max(HPrD_UL, na.rm=TRUE))) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_UL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_line(ggplot2::aes(x=Indices, y=HPrD_LL), color="red", linetype="solid", size=1, na.rm = TRUE) +
ggplot2::geom_ribbon(ggplot2::aes(x=Indices, ymin=HPrD_UL, ymax=HPrD_LL, fill=TRUE), alpha=0.25, show.legend=FALSE) +
ggplot2::scale_fill_manual(values=c("TRUE"="green")) +
ggplot2::geom_point(ggplot2::aes(group=Indices, shape=as.factor(TypeOfdata), color=as.factor(Alarms), stroke = 1.5), show.legend=FALSE, na.rm = TRUE) +
ggplot2::scale_color_manual(values=c("black", "red", "red"), na.value = "black") +
ggplot2::scale_linetype_manual(values=c("Historical"="dotted", "Current"="solid")) +
ggplot2::scale_shape_manual(values=c("Historical"=1, "Current"=19)) +
ggplot2::coord_cartesian(ylim = AdjustedYlim) +
ggplot2::labs(title = main, x = xlab, y = ylab) +
ggplot2::theme(legend.position = "top",
legend.title = ggplot2::element_blank(),
axis.line = ggplot2::element_line(colour = "black", size = 0.5, linetype = "solid"),
panel.background = ggplot2::element_blank(),
plot.title = ggplot2::element_text(hjust = 0.5)
)
print(PCC_historical)
} else { print(PCC) }
}
# List of results
if ( summary_list ) { print(PCC_summary) }
# List of results return in pdf
if ( pdf_report ) {
# save pdf
Cairo::CairoPDF(
paste0( "PCC_results_", paste0( unlist(strsplit(date(), " "))[c(1,2,3,5)], collapse = "_" ), "_",
paste0( unlist(strsplit( unlist(strsplit(date(), " "))[4], ":" )), collapse = "." ),
".pdf" ),
height = 8.264, width = 11.694)
# PCC plot on pdf
print(PCC)
# Prior Posterior plot on pdf
print(PrPostPLOT)
# Results matrix on pdf
# Chunk of code to split results matrix to different pages - Set a default number based on pdf height/width
NRowsPerPage <- 25
if(NRowsPerPage > nrow(PCC_summary)){ FloatingRow <- nrow(PCC_summary) } else { FloatingRow <- NRowsPerPage }
sapply(1:ceiling(nrow(PCC_summary)/NRowsPerPage), function(index) {
if (index==1) { StartingRow <<- 1 }
grid::grid.newpage()
gridExtra::grid.table(PCC_summary[StartingRow:FloatingRow, ])
StartingRow <<- FloatingRow + 1
if( sum(NRowsPerPage, FloatingRow) < nrow(PCC_summary)){ FloatingRow <<- NRowsPerPage + FloatingRow } else { FloatingRow <<- nrow(PCC_summary) }
})
grDevices::dev.off()
}
}
#out <- rgamma(n = 30, shape = 4, rate = 4)
#out[15] <- out[15] + 1.5
#PCC_Gamma(out, al = 4, PriorPosterior_PLOT = TRUE, PCC_PLOT = FALSE)
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.