This document contains the R code for generating estimates of the Upper and Lower Control limits of the process mean with a user specified type-1 error probability based on a truncated saddlepoint approximaiton with a re-expressed fourth cumulant. The first three cumulants of the process distribution are assumed to be known or accurately estimated from data. This code and the equations reference herein correspond to the paper “Estimating \(\bar{X}\) Statistical Control Limits for any Arbitrary Probability Distribution using Re-Expressed Truncated Cumulants” by Braden, P., Chen, B., Matis, T. and Benneyan, J.

####R Input

Load the following libraries into R

#install.packages("pracma")
#install.packages("rootSolve")
library(pracma)
library(rootSolve)

Specify the known/estimated mean, variance, and skewness of the process variable, the subgroup size \(n\), and the desired type-1 error probability \(\alpha\). In the example below, the mean is 12, variance is 36, skewness is 1, n is 10, and alpha is 0.0027

#user-specified input parameters
mean<-12
variance<-36
skewness<-1
n<-10
alpha<-0.0027

From this input, the first three cumulants of the distribution of the process mean are computed.

k1<-mean
k2<-variance/n
k3<-(skewness*variance^(3/2))/n^2

The functions corresponding to \(\theta_0\) in Eq(13), the cumulant generating function (cgf) in Eq(10), the second derivative of the cgf in Eq(11), and \(w\), \(u\), and the cdf of \(\bar{X}\) in Eq(4) are computed.

theta<-function(k1,k2,k3,x) (nthroot((4*k2*k3^3*(2*k2^2-3*k3*(k1-x))),3)-2*k2*k3)/k3^2

K<-function(x) k1*theta(k1,k2,k3,x)+(k2/2)*theta(k1,k2,k3,x)^2+(k3/6)*theta(k1,k2,k3,x)^3+(k3^2/(48*k2))*theta(k1,k2,k3,x)^4

d2K<-function(x) k2+k3*theta(k1,k2,k3,x)+(k3^2/(4*k2))*theta(k1,k2,k3,x)^2

w<-function(x) sign(theta(k1,k2,k3,x))*sqrt(2*(theta(k1,k2,k3,x)*x-K(x)))
u<-function(x) theta(k1,k2,k3,x)*sqrt(d2K(x))
CDF<-function(x) pnorm(eval(w(x)))+dnorm(eval(w(x)))*(1/u(x)-1/w(x))

The Upper Control Limit (UCL) is found by numerically find the root of the cdf evaluated at \(1-\alpha/2\) using the command uniroot in the rootSolve package, for which a user-specified search interval is required. The left endpoint of the search interval should be to the right of the mean, and the right endpoint should be reasonably specified. In this example, the search interval for the root/UCL is specified as \([15,25]\)

#user-specified search interval for UCL
left<-15
right<-25

CDF_Upper<-function(x) CDF(x)-(1-alpha/2)
cat("UCL is:", uniroot(CDF_Upper,c(left,right))$root)
## UCL is: 18.53245

The Lower Control Limit (LCL) is found by numerically find the root of the cdf evaluated at \(\alpha/2\) using the command uniroot in the rootSolve package, for which a user-specified search interval is required. The left endpoint of the search interval should be to the right of zero, and the right endpoint should be left of the mean. In this example, the search interval for the root/UCL is specified as \([0.01,10]\)

#user-specified search interval for LCL
left<-0.1
right<-10

CDF_Lower<-function(x) CDF(x)-(alpha/2)
cat("LCL is:", uniroot(CDF_Lower,c(left,right))$root)
## LCL is: 7.693383

####Complete Unevaluated Code The complete unevaluted code is given below

#install.packages("pracma")
#install.packages("rootSolve")
library(pracma)
library(rootSolve)

#user-specified input parameters
mean<-12
variance<-36
skewness<-1
n<-10
alpha<-0.0027
#
k1<-mean
k2<-variance/n
k3<-(skewness*variance^(3/2))/n^2
theta<-function(k1,k2,k3,x) (nthroot((4*k2*k3^3*(2*k2^2-3*k3*(k1-x))),3)-2*k2*k3)/k3^2
K<-function(x) k1*theta(k1,k2,k3,x)+(k2/2)*theta(k1,k2,k3,x)^2+(k3/6)*theta(k1,k2,k3,x)^3+(k3^2/(48*k2))*theta(k1,k2,k3,x)^4
d2K<-function(x) k2+k3*theta(k1,k2,k3,x)+(k3^2/(4*k2))*theta(k1,k2,k3,x)^2
w<-function(x) sign(theta(k1,k2,k3,x))*sqrt(2*(theta(k1,k2,k3,x)*x-K(x)))
u<-function(x) theta(k1,k2,k3,x)*sqrt(d2K(x))
CDF<-function(x) pnorm(eval(w(x)))+dnorm(eval(w(x)))*(1/u(x)-1/w(x))

#user-specified search interval for UCL
left<-15
right<-25
#
CDF_Upper<-function(x) CDF(x)-(1-alpha/2)
cat("UCL is:", uniroot(CDF_Upper,c(left,right))$root)

#user-specified search interval for LCL
left<-0.1
right<-10
#
CDF_Lower<-function(x) CDF(x)-(alpha/2)
cat("LCL is:", uniroot(CDF_Lower,c(left,right))$root)