R/one_sided_tolerance_interval_tests.R

In stephanGit/leistungstests: Supplies Functions to calculate 'leistungstests'

# kc ------
#' sd-inflation-factor for one sided tolerance intervals.
#'
#'  \code{kc} calculates the factor for one sided tolerance intervalls
#'  according to the Formula A.13 in ISO 16269-6.
#'
#' @aliases kc ltl
#'
#' @param p the desired coverage probability of the tolerance intervall
#' @param alpha Niveau of the test, i.e. Confidence-level = 1 - alpha
#' @param n the sample size
#'
#' @return the factor kc in x-kc*s or x+kc*x
#'
#' @seealso \url{https://www.iso.org/obp/ui/#iso:std:iso:16269:-6:ed-2:v1:en} for an intro
#' to the iso-norm.
#'
#' @export
kc <- function(p, alpha, n) {
  stats::qt(p = (1 - alpha), df = (n - 1), ncp = sqrt(n) * stats::qnorm(p = p))/base::sqrt(n)
}


# p_reject_ostit --------------------------
#' Probability to reject one-sided-tolerance-interval test
#'
#' @description
#' Calculates the rejection-probability for one-sided-tolerance-interval tests
#' given the assumed process-mean and process-sd for specified test-parameters
#'
#' @param process_mean assumed process-mean
#' @param process_sd assumed process-sd
#' @param sample_size sample-size
#' @param spec_limit specification limit
#' @param spec_type ("lsl" or "usl") specification-type, i.e. lower or upper one-sided spec-limit
#' @param p_min desired minimal coverage-probability of tolerance-intervall
#' @param power desired test-power
#'
#' @details
#' A one-sided-tolerance interval test is conducted by calculating a tolerance-intervall
#' according to ISO 16269-6. The test is passed if this tolerance-intervall is within the
#' specification limits. The desired power of the test is exactly the confidence-level for
#' the tolerance-interval!
#'
#' @return Probabilty to reject H_0 given the process-mean, standard-deviation and specification
#'
#' @export
p_reject_ostit <- function(process_mean,
                           process_sd,
                           sample_size,
                           spec_limit,
                           spec_type,
                           p_min,
                           power){

  slope_factor <- ((process_mean - spec_limit) / process_sd)
  if ( spec_type == "usl") {
    slope_factor <- -1 * slope_factor
  }
  p_reject <- pt(q = sqrt(sample_size) * kc(p = p_min, alpha = 1-power, n = sample_size),
                 df = (sample_size - 1),
                 ncp = sqrt(sample_size) * slope_factor)

  return(p_reject)
}


# sample_size_hint_ostit --------------------
#' Suggest sample-size for one-sided tolerance intervall test
#'
#' @description
#' Calculates a sample-size to achieve a desired probabilty to pass the ostit
#' for specified lower tolerance-limit and (hypothized) process-mean and process-sd
#' and tolerance-intervall parameters coverage-probability and alpha.
#'
#' @details
#' sample-size wird bestimmt indem die Nullstelle der Funktion:
#' p_reject_ostit(sample_size) - (1-p_success) für sample_size zwischen 2 und 1000 gesucht wird.
#' Ist der Prozess so gut das p_reject ~ 0 für alle n=2..1000 so wird sample_size = 3 zurückgegeben
#'
#' Ist der Prozess i.o. aber dicht an der Grenze, so kann p_success maximal (1-power) sein.
#'
#' @param p_success desired minimal probability to pass the ostit
#' @param process_mean assumed process mean (estimated from Vorentwicklungs-Versuchen)
#' @param process_sd assumed process sd (estimated from Vorentwicklungs-Versuchen)
#' @param spec_limit spec limit for variable
#' @param spec_type ("lsl" or "usl") specification-type, i.e. lower or upper one-sided spec-limit
#' @param p_min geforderte minimale i.o. rate des prozesses
#' @param power geforderte Power des Tests
#'
#' @return sample-size
#'
#' @export
sample_size_hint_ostit <- function(p_success,
                                   process_mean,
                                   process_sd,
                                   spec_limit,
                                   spec_type,
                                   p_min,
                                   power){

  if (is.na(p_success) | is.na(process_mean) | is.na(process_sd) | is.na(spec_limit) |
      is.na(spec_type) | is.na(p_min) | is.na(power)){
    warning("At least on arg is NA which must not be!")
    return(NA)
  }

  if (process_sd == 0){
    warning(("Process-SD = 0; return NA"))
    return(NA)
  }

  if (p_success < 1 - power){
    warning("Erfolgswahrscheinlichkeit kann nicht kleiner sein als 1-power")
    return(NA)
  }

  if (process_mean < spec_limit + qnorm(p_min) * process_sd & spec_type == "lsl"){
    warning("Process erfuellt die Anforderungen nicht!")
    return(NA)
  }

  if (process_mean > spec_limit - qnorm(p_min) * process_sd & spec_type == "usl"){
    warning("Process erfuellt die Anforderungen nicht!")
    return(NA)
  }

  score_function <- function(sample_size){
    p_reject_ostit(
      process_mean = process_mean,
      process_sd = process_sd,
      spec_limit = spec_limit,
      spec_type = spec_type,
      p_min = p_min,
      sample_size = sample_size,
      power = power
    ) - (1 - p_success)
  }

  sample_size_limits <- c(2, 1000)

  if ( score_function(sample_size_limits[1]) < 0 &
       score_function(sample_size_limits[2]) < 0 &
       (
         (process_mean > spec_limit & spec_type == "lsl") |
        (process_mean < spec_limit & spec_type == "usl")
       )
  ){
    return(3)
  }

  if ( score_function(sample_size_limits[1]) > 0 &
       score_function(sample_size_limits[2]) > 0 ){
    warning(paste("Keine Stichprobe kleiner ", sample_size_limits[2],
                  " kann die Anforderung erfuellen! Rueckgabewert auf 10000 gesetzt! ", sep = ""))
    return(10000)
  }

  return(ceiling(uniroot(score_function, interval = c(sample_size_limits[1],sample_size_limits[2]))$root))
}