R/mewma.arl.f.R

In spc: Statistical Process Control -- Calculation of ARL and Other Control Chart Performance Measures

# Computation of MEWMA ARLs (multivariate mean monitoring), returns function
mewma.arl.f <- function(l, cE, p, delta=0, r=20, ntype=NULL, qm0=20, qm1=qm0) {
  if ( l<=0 | l>1 )		stop("l has to be between 0 and 1")
  if ( cE<=0 )			stop("threshold c has to be positive")
  if ( p<1 )			stop("wrong dimension parameter")
  if ( delta<0 )		stop("wrong magnitude value")
  
  if ( r<4 )			stop("resolution too small")
  if ( qm0<5 )			stop("more quadrature nodes needed")
  if ( qm1<5 )			stop("more quadrature nodes needed")
  
  if ( is.null(ntype) ) {
    if ( delta <1e-10 ) {
      ntype <- "gl2"
    } else {
      #if ( p %in% c(2,4) ) {
      if ( p==2 ) {
        ntype <- "gl3"
      } else {
        ntype <- "gl5"
      }
    }
  }
  
  # collocation basis of Chebshev polynomials
  Tn <- Vectorize(function(z, n) {
    if ( n==0 ) result <- 1
    if ( n==1 ) result <- z
    if ( n==2 ) result <- 2*z^2 - 1
    if ( n==3 ) result <- 4*z^3 - 3*z
    if ( n==4 ) result <- 8*z^4 - 8*z^2 + 1
    if ( n==5 ) result <- 16*z^5 - 20*z^3 + 5*z
    if ( n>5 )  result <- cos( n*acos(z) ) 
    result
  })
  
  qtyp <- pmatch(tolower(ntype), c("gl", "co", "ra", "cc", "mc", "sr", "co2", "gl2", "gl3", "gl4", "gl5", "co3", "co4", "ngl1", "ngl2", "ngl3", "ngl4", "ngl5")) - 1
  if ( is.na(qtyp) )		stop("invalid type of numerical algorithm")

  if ( abs(delta) < 1e-10 ) { # in-control    
    LENGTH <- 3*r
    zeug <- .C("mewma_arl_f", as.double(l), as.double(cE), as.integer(p), as.double(delta),
                              as.integer(r), as.integer(qtyp), as.integer(qm0), as.integer(qm1),
                              ans=double(length=LENGTH), PACKAGE="spc")$ans
                              
    g <- zeug[1:r]                          
    w <- zeug[1:r + r]
    z <- zeug[1:r + 2*r]
    
    # helper functions
    cE <- cE * l/(2-l)
    l2 <- ( (1-l)/l )^2 
    fchi <- function(a, u) dchisq( u/l^2, p, ncp=l2*a ) / l^2
    FCHI <- function(a, u) pchisq( u/l^2, p, ncp=l2*a )
    
    if ( qtyp %in% c(0,2,5) ) arl <- Vectorize(function(a) 1 + sum( w * fchi(a, z) * g ), "a") # ordinary GL or Radau or Simpson rule Nystroem 
    if ( qtyp==7 ) arl <- Vectorize(function(a) 1 + sum( w * fchi(a, z^2) * g * 2*z ), "a") # GL Nystroem with ()^2 substitution
    if ( qtyp==1 ) arl <- Vectorize(function(a) sum( Tn( (2*a-cE)/cE, 0:(r-1) ) * g ), "a") # collocation
    if ( qtyp==3 ) arl <- Vectorize(function(a) 1 + sum( w * fchi(a, z^2) * g ) * cE/2 , "a") # Clenshaw-Curtis
    if ( qtyp==4 ) arl <- Vectorize(function(a) 1 + sum( ( FCHI(a, z^2) - FCHI(a, c(0, z[-length(z)]^2)) ) * g), "a") # Markov chain (Runger/Prabhu)
    
  } else { # out-of-control
    if ( qtyp==4 ) {
      cE_ <- sqrt( cE * l/(2-l) )
      w  <- 2*cE_/( 2*r + 1 )
      ii <- function(ix, iy) (ix-r)^2 + iy^2 < cE_^2/w^2
      CIRC <- as.vector( t(outer( 0:(2*r), 0:r, ii)) )
      dQ <- sum(CIRC)
      LENGTH <- dQ
    } else {
      r2 <- r^2
      LENGTH <- r2 + 4*r
    }
    
    zeug <- .C("mewma_arl_f", as.double(l), as.double(cE), as.integer(p), as.double(delta),
                              as.integer(r), as.integer(qtyp), as.integer(qm0), as.integer(qm1),
                              ans=double(length=LENGTH), PACKAGE="spc")$ans                              
                              
    if ( qtyp!=4 ) {
      g  <- zeug[1:r2]
      w0 <- zeug[1:r + r2]
      z0 <- zeug[1:r + r2 + r]
      w1 <- zeug[1:r + r2 + 2*r]
      z1 <- zeug[1:r + r2 + 3*r]
    } else {
      g <- zeug[1:dQ]
    }
    
    # helpers
    l2 <- ( (1-l)/l )^2
    h <- cE * l/(2-l)
    rdc <- l * sqrt( delta/h )
    sig <- l / sqrt( h )
    
    lsd  <- l * sqrt( delta )    
    
    if ( qtyp %in% c(0, 2, 3, 5) ) arl <- Vectorize(function(a, b) { # ordinary GL or Radau or CC or Simpson rule Nystroem 
      if ( abs(h-a) < 1e-10 ) a_ <- 1 else a_ <- ( a - b^2/delta ) / ( h - b^2/delta )
      b_ <- b / sqrt( delta * h )
      m <- rdc + (1-l)*b_
      eta <- l2 * h * (1-b_^2) * a_
      if ( eta < 1e-10 ) eta <- 0
      result <- 1
      for ( i in 1:r ) {
        korr <- h * (1-z1[i]^2) / l^2 
        outer <- korr * w1[i] * dnorm( z1[i], mean=m, sd=sig)
        inner <- sum( w0 * dchisq( korr*z0, p-1, eta) * g[ (i-1)*r + 1:r ] )
        result <- result + inner * outer
      }
      result
    })
    
    if ( qtyp==7 ) arl <- Vectorize(function(a, b) { # GL Nystroem with ()^2 substitution
      if ( abs(h-a) < 1e-10 ) a_ <- 1 else a_ <- ( a - b^2/delta ) / ( h - b^2/delta )
      b_ <- b / sqrt( delta * h )
      m <- rdc + (1-l)*b_
      eta <- l2 * h * (1-b_^2) * a_
      if ( eta < 1e-10 ) eta <- 0
      result <- 1
      for ( i in 1:r ) {
        korr <- h * (1-z1[i]^2) / l^2 
        outer <- korr * w1[i] * dnorm( z1[i], mean=m, sd=sig)
        inner <- sum( w0 * dchisq( korr*z0^2, p-1, eta) * 2*z0 * g[ (i-1)*r + 1:r ] )
        result <- result + inner * outer
      }
      result
    })
        
    if ( qtyp==8 ) arl <- Vectorize(function(a, b) { # GL Nystroem with ()^2 plus sin() substitution
      if ( abs(h-a) < 1e-10 ) a_ <- 1 else a_ <- ( a - b^2/delta ) / ( h - b^2/delta )
      b_ <- b / sqrt( delta * h )
      m <- rdc + (1-l)*b_
      eta <- l2 * h * (1-b_^2) * a_
      if ( eta < 1e-10 ) eta <- 0
      result <- 1
      for ( i in 1:r ) {
        korr <- h * ( 1 - sin(z1[i])^2 ) / l^2 
        outer <- korr * w1[i] * dnorm( sin(z1[i]), mean=m, sd=sig) * cos(z1[i])
        inner <- sum( w0 * dchisq( korr*z0^2, p-1, eta) * 2*z0 * g[ (i-1)*r + 1:r ] )
        result <- result + inner * outer
      }
      result
    })
    
    if ( qtyp==9 ) arl <- Vectorize(function(a, b) { # GL Nystroem with ()^2 plus tan() substitution
      if ( abs(h-a) < 1e-10 ) a_ <- 1 else a_ <- ( a - b^2/delta ) / ( h - b^2/delta )
      b_ <- b / sqrt( delta * h )
      m <- rdc + (1-l)*b_
      eta <- l2 * h * (1-b_^2) * a_
      if ( eta < 1e-10 ) eta <- 0
      result <- 1
      for ( i in 1:r ) {
        korr <- h * ( 1 - tan(z1[i])^2 ) / l^2 
        outer <- korr * w1[i] * dnorm( tan(z1[i]), mean=m, sd=sig) / cos(z1[i])^2
        inner <- sum( w0 * dchisq( korr*z0^2, p-1, eta) * 2*z0 * g[ (i-1)*r + 1:r ] )
        result <- result + inner * outer
      }
      result
    })   
 
    if ( qtyp==10 ) arl <- Vectorize(function(a, b) { # GL Nystroem with ()^2 plus sinh() substitution
      norm <- sinh(1)
      if ( abs(h-a) < 1e-10 ) a_ <- 1 else a_ <- ( a - b^2/delta ) / ( h - b^2/delta )
      b_ <- b / sqrt( delta * h )
      m <- rdc + (1-l)*b_
      eta <- l2 * h * (1-b_^2) * a_
      if ( eta < 1e-10 ) eta <- 0
      result <- 1
      for ( i in 1:r ) {
        korr <- h * ( 1 - (sinh(z1[i])/norm)^2 ) / l^2         
        outer <- korr * w1[i] * dnorm( sinh(z1[i])/norm, mean=m, sd=sig) * cosh(z1[i])/norm
        inner <- sum( w0 * dchisq( korr*z0^2, p-1, eta) * 2*z0 * g[ (i-1)*r + 1:r ] )
        result <- result + inner * outer
      }
      result
    })
    
    if ( qtyp %in% c(1, 6, 11, 12) ) arl <- Vectorize(function(a, b) { # collocation 
      if ( abs(h-a) < 1e-10 ) a_ <- 1 else a_ <- ( a - b^2/delta ) / ( h - b^2/delta )
      b_ <- b / sqrt( delta * h )
      result <- 0
      for ( i in 1:r ) {
        outer <- Tn( 2*a_-1, i-1 ) 
        inner <- sum( Tn( b_, 0:(r-1)) * g[ (i-1)*r + 1:r ] )
        result <- result + inner * outer
      }
      result
    })
    
    if ( qtyp==4 ) arl <- Vectorize(function(a, b) { # Markov chain (Runger/Prabhu)
      #a   <- sqrt(a)
      cE_ <- sqrt( cE * l/(2-l) )
      w  <- 2*cE_/( 2*r + 1 )
      wl <- w^2 / l^2
      ii <- function(ix, iy) (ix-r)^2 + iy^2 < cE_^2/w^2
      Vf <- function(iy,jy) pchisq( (jy+.5)^2*wl, p-1, ncp=(iy*w)^2*l2) - as.numeric(jy>0)*pchisq( (jy-.5)^2*wl, p-1, ncp=(iy*w)^2*l2)
      Hf <- function(ix,jx) pnorm( (-cE_+(jx+1)*w-(1-l)*(-cE_+(ix+.5)*w) )/l, mean=delta) - pnorm( (-cE_+jx*w-(1-l)*(-cE_+(ix+.5)*w) )/l, mean=delta)
      CIRC <- as.vector( t(outer( 0:(2*r), 0:r, ii)) )
      Vv <- Vf( a/w, 0:r )
      Hv <- Hf( (b+cE_)/w-.5, 0:(2*r) )
      dQ <- sum(CIRC)
      Qv <- as.vector( Vv %o% Hv )
      Qv_ <- Qv[ CIRC ]
      result <- 1 + sum( Qv_ * g )
    })
    
    if ( qtyp %in% c(13, 14, 15, 16, 17) ) arl <- Vectorize(function(a, b) { # new GL designs
      gam <- 0
      #if ( a > 0 ) gam <- b / sqrt( a * delta )
      if ( a > 0 ) gam <- b / sqrt( a )
      mij <- lsd + (1 - l ) * sqrt(a) * gam
      ncpij <- l2 * a * ( 1 - gam^2 )
      norm <- sinh(1)
      result <- 1
      for ( i in 1:r ) {
        if ( qtyp == 13 ) korr <- w0[i] * sqrt( z0[i] ) / l^2
        if ( qtyp > 13 ) korr <- 2 * w0[i] * z0[i]^2 / l^2
        wl <- z1
        korr2 <- rep(1, r)
        if ( qtyp == 15 ) { wl <- sin( z1 ); korr2 <- cos( z1 ) }
        if ( qtyp == 16 ) { wl <- tan( z1 ); korr2 <- 1 / ( cos( z1 )^2 ) }
        if ( qtyp == 17 ) { wl <- sinh( z1 ) / norm;  korr2 <- cosh( z1 ) / norm }
        if ( qtyp == 13 ) term <- korr * sum( w1 * dnorm( (sqrt(z0[i])*wl - mij)/l ) / l * dchisq( z0[i]*( 1 - wl^2 ) / l^2, p-1, ncpij ) * g[ (i-1)*r + 1:r ] * korr2 )
        if ( qtyp > 13 )  term <- korr * sum( w1 * dnorm( (z0[i]*wl - mij)/l ) / l * dchisq( z0[i]^2*( 1 - wl^2 ) / l^2, p-1, ncpij ) * g[ (i-1)*r + 1:r ] * korr2 )
        result <- result + term 
      }
      result
    })
  } 
  
  arl
}