[R] qcc
Leif Kirschenbaum
leif at reflectivity.com
Wed Nov 30 19:37:41 CET 2005
If you examine the code for the function "violating.runs" in the qcc package (try "violating.runs") you can deconstruct it to find that the code flags runs where there are .qcc.options$run.length or more points in a row on one side of the center (process mean).
However, the classic Shewhart rules dictate that a run of monotonically increasing or decreasing points of <run length> is also a run violation. I include here my modified version of the function "violating.runs" which also checks for these violations:
##
## Correct some typos in violating.runs from qcc package
## Added test for run.length of points monotonically increasing or decreasing
## The simplest way is to re-run the code but with "diffs"
## representing the sign of the difference from one point to the next
##
violating.runs<-function (object, run.length = qcc.options("run.length"))
{
center <- object$center
statistics <- c(object$statistics, object$new.statistics)
cl <- object$limits
violators <- numeric()
for(i in 1:2){
diffs <- statistics - center
if(i==2) {
diffs <- c(0,diff(statistics))
## need to decrement run.length since we're looking at differences between points
run.length<-run.length-1
}
diffs[diffs > 0] <- 1
diffs[diffs < 0] <- -1
runs <- rle(diffs)
index.lengths <- (1:length(runs$lengths))[runs$lengths >= run.length]
index.stats <- 1:length(statistics)
vruns <- rep(runs$lengths >= run.length, runs$lengths)
vruns.above <- (vruns & (diffs > 0))
vruns.below <- (vruns & (diffs < 0))
rvruns.above <- rle(vruns.above)
rvruns.below <- rle(vruns.below)
vbeg.above <- cumsum(rvruns.above$lengths)[rvruns.above$values]-(rvruns.above$lengths - run.length)[rvruns.above$values]
vend.above <- cumsum(rvruns.above$lengths)[rvruns.above$values]
vbeg.below <- cumsum(rvruns.below$lengths)[rvruns.below$values]-(rvruns.below$lengths - run.length)[rvruns.below$values]
vend.below <- cumsum(rvruns.below$lengths)[rvruns.below$values]
if (length(vbeg.above)) {
for (i in 1:length(vbeg.above)) violators <- c(violators, vbeg.above[i]:vend.above[i])
}
if (length(vbeg.below)) {
for (i in 1:length(vbeg.below)) violators <- c(violators, vbeg.below[i]:vend.below[i])
}
} ## ENDOF for i in 1:2
return(violators)
}
> 3)Is there any more criterias made somewhere ?
The other criteria is found by the function "beyond.limits" which checks to see if any points are beyond the upper or lower control limits.
I believe that Prof. Scrucca wrote the qcc package as a tool to demonstrate process control to his students taking his statistics classes. Hence the statistics are excellent. For example note the use of log-gamma and exponential functions in the function "sd.xbar" to calculate the constant "c4", which avoids the divison of extremely large [~ 1E20] numbers:
c4 <- function(n) sqrt(2/(n-1)) * exp(lgamma(n/2) - lgamma((n-1)/2)))
where a textbook might write it as:
c4 <- function(n) sqrt(2/(n-1)) * (gamma(n/2) / gamma((n-1)/2)))
Thus I suggest that although qcc provides an excellent basis for constructing statistical control reports, that you should be prepared to modify it for your purposes (as I am heavily doing).
P.S. I change .qcc.options using syntax such as ".qcc.options$violating.runs<-7".
Leif S. Kirschenbaum, Ph.D.
Senior Yield Engineer
Reflectivity, Inc.
408-737-8100 x307
408-737-8153 (Fax)
> Date: Tue, 29 Nov 2005 02:08:57 +0200
> From: Tommi Viitanen <totavi at utu.fi>
> Subject: [R] qcc
>
> violating.runs
[deleted]
>
> that the criteria for the violating is 5 but
> 1)I cannot find "5" in the code of the function. Where is the "5" ?
> 2)What is the easiest way to change it ?
> 3)Is there any more criterias made somewhere ?
>
> Yours sincerelly, Tommi Viitanen
>
[deleted]
> ------------------------------
>
> Date: Tue, 29 Nov 2005 08:58:24 +0100
> From: Uwe Ligges <ligges at statistik.uni-dortmund.de>
> Subject: Re: [R] qcc
> To: Tommi Viitanen <totavi at utu.fi>
>
>
> See ?violating.runs:
> violating.runs(object, run.length = qcc.options("run.length"))
> ^^^^^ ^^^ ^^^ ^^ ^^^ ^^^
>
[deleted]
>
> Uwe Ligges
More information about the R-help
mailing list