# R/code_generic.R In andreaskiermeier/AcceptanceSampling: Creation and Evaluation of Acceptance Sampling Plans

#### Documented in find.plan

```## code_generic.R ---
##
## Author: Andreas Kiermeier
##
## Created: 21 Aug 2007
##
## Purpose: Generic code (mostly) which applies to all types of sampling plans
##
## Changes:
## 21Aug07: * Created
## 29Dec16: * find.plan: Included check for missing N for hypergeometric
## 06Apr22: * find.k: Include interval=c(0,1000) in call to find.k
## ----------------------------------------------------------------------

setGeneric("assess", function(object, PRP, CRP, print=TRUE)
standardGeneric("assess"))

check.paccept <-
function(pa){
## Purpose: Utility function to check that supplied P(accept) values
##          fall within [0,1]
## ----------------------------------------------------------------------
## Arguments:
## pa: a vector of P(accept) values
## ----------------------------------------------------------------------
## Author: Andreas Kiermeier, Date: 16 May 2007, 10:19
if (any(pa < 0) | any(pa > 1))
return(FALSE)
return(TRUE)
}

check.quality <-
function(pd, type){
## Purpose: Utility function to check that supplied Proportion defective
##          values fall within
##          [0,1] for the binomial or hypergeometric
##          [0,inf] for the poisson
## ----------------------------------------------------------------------
## Arguments:
## pd: a vector of proportion defective values
## ----------------------------------------------------------------------
## Author: Andreas Kiermeier, Date: 16 May 2007, 10:19
if (any(pd < 0))
return(FALSE)
if (type %in% c("binomial", "hypergeom", "normal") & any(pd > 1))
return(FALSE)
return(TRUE)
}

## Utility to find k for a given n
find.k <- function(n, pd, pa, interval=c(0,5)){
tmp <- uniroot(function(x, n, pd, pa){
pt(x*sqrt(n), df=n-1, ncp=-qnorm(pd)*sqrt(n)) - pa},
interval=interval, n=n, pd=pd, pa=1-pa)
return(tmp\$root)
}

find.plan <- function(PRP, CRP,
type=c("binomial","hypergeom","poisson","normal"),
N,
s.type=c("known", "unknown"))
{
## Purpose: Find the sampling plan with the smallest sample size, which
##          meets a prespecified Producer and Consumer Risk Points.
##
##          The convention used here, as in many books, is to use equality
##          for the Producer Risk Point rather than the consumer risk point.
##
##          No consideration is given to "cost functions".
## ----------------------------------------------------------------------
## Arguments:
## PRP   : Producer risk point in the form c(pdefect, paccept)
## CRP   : Consumer risk point in the form c(pdefect, paccept)
## N     : Population size - only used for hypergeomtric distribution
## type  : The distributional assumption
## s.type: Only used for 'normal' distribution - indicates whether the
##         standard deviation is known or unknown (use sample s.d.)
## ----------------------------------------------------------------------
## Author: Andreas Kiermeier, Date: 20 Aug 2007, 12:09

type <- match.arg(type)
s.type <- match.arg(s.type)

## Check that N is supplied for hypergeometric distribution
if(type=="hypergeom" & missing(N))
stop("N must be supplied for the hypergeometric distribution.")
## Needs checking that risk points are "valid" - use existing functions
if (missing(PRP) | missing(CRP))
stop("Poducer and Consumer Risk Points must be provided.")
else if(!check.quality(PRP[1], type=type) |
!check.paccept(PRP[2]) )
stop("Producer Risk Point - Quality and/or desired P(accept) out of bounds")
else if(!check.quality(CRP[1], type=type) |
!check.paccept(CRP[2]) )
stop("Consumer Risk Point - Quality and/or desired P(accept) out of bounds")
else if(CRP[1] <= PRP[1])
stop("Consumer Risk Point quality must be greater than Producer Risk Point quality")

## Attributes Sampling Plan - Binomial distribution
if (type == "binomial") {
c <- 0
n <- c+1
repeat {
if (calc.OCbinomial(n=n,c=c,r=c+1,pd=CRP[1]) > CRP[2])
n <- n + 1
else if (calc.OCbinomial(n=n,c=c,r=c+1,pd=PRP[1]) < PRP[2])
c <- c + 1
else
break
}
return(list(n=n, c=c, r=c+1))
}
## Attributes Sampling Plan - Hypergeometric distribution
if (type == "hypergeom") {
c <- 0
n <- c+1
repeat {
if (calc.OChypergeom(n=n,c=c,r=c+1,N=N,D=CRP[1]*N) > CRP[2])
n <- n + 1
else if (calc.OChypergeom(n=n,c=c,r=c+1,N=N,D=PRP[1]*N) < PRP[2])
c <- c + 1
else
break
}
return(list(n=n, c=c, r=c+1))
}
## Attributes Sampling Plan - Poisson distribution
if (type == "poisson") {
c <- 0
n <- c+1
repeat {
if (calc.OCpoisson(n=n,c=c,r=c+1,pd=CRP[1]) > CRP[2])
n <- n + 1
else if (calc.OCpoisson(n=n,c=c,r=c+1,pd=PRP[1]) < PRP[2])
c <- c + 1
else
break
}
return(list(n=n, c=c, r=c+1))
}
## Variables Sampling Plan - Normal distribution
else if (type=="normal") {
## With known standard deviation
if (s.type=="known") {
n <- ceiling( ((qnorm(1-PRP[2]) + qnorm(CRP[2]))/
(qnorm(CRP[1])-qnorm(PRP[1])) )^2)
k <- qnorm(1-PRP[2])/sqrt(n) - qnorm(PRP[1])
return(list(n=n, k=k, s.type=s.type))
}
## With unknown standard deviation
else if (s.type=="unknown") {
n <- 2 ## Need a minimum of 1 degree of freedom (=n-1) for the NC t-dist
k <- find.k(n, PRP[1], PRP[2], interval=c(0,1000))
pa <- 1- pt(k*sqrt(n), df=n-1, ncp=-qnorm(CRP[1])*sqrt(n))
while(pa > CRP[2]){
n <- n+1
k <- find.k(n, PRP[1], PRP[2], interval=c(0,1000))
pa <- 1-pt(k*sqrt(n), df=n-1, ncp=-qnorm(CRP[1])*sqrt(n))
}
return(list(n=n, k=k, s.type=s.type))
}
}
}

## x1 <- find.plan(c(0.05, 0.95), c(0.15, 0.075), type="bin")
## x <- OC2c(x1\$n, x1\$c, x1\$r, type="bin")
## assess(x, c(0.05, 0.95), c(0.15, 0.075))

## x1 <- find.plan(c(0.05, 0.95), c(0.15, 0.075), type="hyp", N=100)
## x <- OC2c(x1\$n, x1\$c, x1\$r, type="hyp", N=100)
## assess(x, c(0.05, 0.95), c(0.15, 0.075))

## x1 <- find.plan(c(0.05, 0.95), c(0.15, 0.075), type="pois")
## x <- OC2c(x1\$n, x1\$c, x1\$r, type="pois")
## assess(x, c(0.05, 0.95), c(0.15, 0.075))

## The following examples come from Guenther's book

## PRP <- c(0.01, 0.95)
## CRP <- c(0.10, 0.1)
## x1 <- find.plan(PRP=PRP, CRP=CRP, type="nor", s.type="unknown")
## x <- OCvar(x1\$n, x1\$k, s.type=x1\$s.type, pd=seq(0,0.2, by=0.002))
## plot(x)
## points(PRP[1], PRP[2], col="red", pch=19); points(CRP[1], CRP[2], col="red", pch=19)
## assess(x, PRP=PRP, CRP=CRP)

## PRP <- c(0.05, 0.95)
## CRP <- c(0.20, 0.1)
## x1 <- find.plan(PRP=PRP, CRP=CRP, type="nor", s.type="known")
## x <- OCvar(x1\$n, x1\$k, s.type=x1\$s.type, pd=seq(0,0.2, by=0.002))
## plot(x)
## points(PRP[1], PRP[2], col="red", pch=19); points(CRP[1], CRP[2], col="red", pch=19)
## assess(x, PRP=PRP, CRP=CRP)

### Local Variables:
### comment-start: "## "
### fill-column: 80
### End:
```
andreaskiermeier/AcceptanceSampling documentation built on April 7, 2022, 12:37 a.m.