Simulaciones

Author

Ana María Valencia

Librerías

library(fda)
library(MFPCA)
library(funData)

CASO UNIVARIADO

En este caso se realizan varias simulaciones de observaciones funcionales para evaluar el desempeño y comportamiento del índice de capacidad multivariado basado en componentes principales dados los límites de especificación preestablecidos. En cada caso, se simulan \(N=250\) observaciones funcionales univariadas tal que:

\[f_i(t) = \sum_{m=1}^M \varepsilon_{im}\phi_m(t)\]

con \(t \in T\) y \(i=1,...,250\).

Adicionalmente, los scores \(\varepsilon_{im}\) son muestras independientes de una distribución \(N(0, \lambda_m)\) para valores propios con decrecimiento lineal, exponencial o correspondientes a un proceso de Wiener.

Se consideran distintas configuraciones para cada escenario, donde varían:

  • Los valores de M
  • Tipo de funciones con que serán simuladas las curvas (“Poly”: Polinomios de Legendre o “Fourier”: Funciones de la base de Fourier).
  • El tipo de decrecimiento de los valores propios: (“linear”: \(\nu_m = \frac{M+1-m}{m}\), “exponential” = \(\nu_m = exp(-\frac{m-1}{2})\), “wiener” = \(\nu_m = \frac{1}{(\pi/2 \cdot (2m-1))^2}\)).
  • El tipo de base con que son ajustadas las observaciones simuladas (Fourier o B-Splines).

Límites de especificación

Se crea un conjunto de datos referencia para establecer los límites de especificación que serán utilizados en los escenarios propuestos de simulación. Para este caso, se consideran \(M=20\) funciones de la base de Fourier en el intervalo \([0,1]\) para el proceso univariado. Los límites de control están dados por:

\[USL(t)=\bar{x}(t)+3.5 \sigma(t)\] \[LSL(t)=\bar{x}(t)-3.5\sigma(t)\]

Code 
# Construcción límites de especificación

set.seed(2143)

t <- seq(0, 1, length.out = 100)

sim_uni <- simFunData(argvals = t, M=20, eFunType = "Fourier", eValType = "exponential", N=250)

mu <- colMeans(sim_uni$simData@X)
sigma <- apply(sim_uni$simData@X, 2, sd)
k <- 3.5

# Límites de especificación
USL_uni <- mu + k * sigma
LSL_uni <- mu - k * sigma

# Gráfico
matplot(t, t(sim_uni$simData@X), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "", ylim = c(-6, 6))
lines(t, USL_uni, col = "red", lwd = 3)
lines(t, LSL_uni, col = "red", lwd = 3)
lines(t, colMeans(sim_uni$simData@X), lwd = 3, col = "blue")

Code 
USL_fun_uni <- funData(argvals = list(t), X = matrix(USL_uni, nrow = 1))
LSL_fun_uni <- funData(argvals = list(t), X = matrix(LSL_uni, nrow = 1))
mean_fun_uni <- funData(argvals = list(t), X = matrix(mu, nrow = 1))
Code 
# ESCENARIO 1: M=3, eFunType = "Fourier", eValType = "exponential"

set.seed(2143)

escenario1 <- simFunData(argvals = t, M=3, eFunType = "Fourier", eValType = "exponential", N=250)
base <- create.fourier.basis(rangeval = c(0,1), nbasis = 20)
fd_obj <- funData2fd(escenario1$simData, basisobj = base)

# PCA
pca <- pca.fd(fd_obj, nharm = 3)
scores <- pca$scores
lambda <- pca$values
ncomp <- 3
T2_uni <- rowSums(sweep(scores[,1:ncomp]^2, 2, lambda[1:ncomp], "/"))

alpha = 1-sqrt(1-0.01)

# Límite de control

UCL_T2_FPCA_sim <- function(simData, alpha = 0.01, ncomp = 3, R = 2000, nbasis = 20){
  n <- nObs(simData)
  t_grid <- simData@argvals[[1]]
  basis_obj <- create.fourier.basis(rangeval = range(t_grid), nbasis = nbasis)
  T2_vmax <- numeric(R)
  for(i in 1:R){
    indices <- sample(1:n, n, replace = TRUE)
    X_boot <- simData@X[indices, ]
    fun_boot <- funData(argvals = list(t_grid), X = X_boot)
    fd_boot <- funData2fd(fun_boot,  basisobj = basis_obj)
    pca_boot <- pca.fd(fd_boot, nharm = ncomp, centerfns = TRUE)
    scores <- pca_boot$scores[,1:ncomp]
    lambda <- pca_boot$values[1:ncomp]
    T2_vals <- rowSums(sweep(scores^2, 2, lambda, "/"))
    T2_vmax[i] <- max(T2_vals)
  }
  quantile(T2_vmax, probs = 1 - alpha)
}

UCL_T2_uni <- UCL_T2_FPCA_sim(simData = escenario1$simData, alpha = alpha, ncomp = 3, R = 3000)

# Carta T2 

par(mfrow=c(1,1))
plot(T2_uni, type = "p", pch = 16, cex = 0.8, xlab = "Observación",
     ylab = expression(T^2), main = expression("Carta " * T^2),
     ylim=c(0,20))
abline(h = UCL_T2_uni, col = "red", lwd = 2, lty = 2)

Code 
# Proceso bajo control


# Carta Q 

t_grid <- sim_uni$simData@argvals[[1]]
dt <- diff(t_grid)[1]

scores <- pca$scores

###### CURVAS GORRO 

gorros <- vector("list", nrow(scores))

coef_harm <- pca$harmonics$coefs   
basis <- pca$harmonics$basis
mean_coef <- pca$meanfd$coefs

n <- nrow(scores)
K <- ncol(scores)

for(i in 1:n){
  coef_i <- rep(0, nrow(coef_harm))
  for(k in 1:K){
    coef_i <- coef_i + scores[i,k] * coef_harm[,k]
  }
  coef_i <- coef_i + mean_coef
  gorros[[i]] <- fd(coef_i, basis)
}

coef_matrix <- do.call(cbind, lapply(gorros, function(f) f$coefs))
gorros_fd <- fd(coef_matrix, basis)

# Reconstrucciones KL

gorros_mat <- sapply(gorros, function(f) eval.fd(t, f))
matplot(t, gorros_mat, type = "l", lty = 1, col = "gray",
        main = "Curvas reconstruidas", xlab = "t", ylab = "X(t)")

# Media funcional
X <- eval.fd(t, fd_obj)
mu_X <- rowMeans(X)

lines(t, mu_X, lwd = 3, col = "black")

Code 
# Estadístico Q
residuos <- fd_obj - gorros_fd
Q <- diag(inprod(residuos, residuos))

UCL_Q_FPCA_sim <- function(simData, alpha = 0.01, ncomp = 3, R = 2000, nbasis = 20){
  n <- nObs(simData)
  t_grid <- simData@argvals[[1]]
  basis_obj <- create.fourier.basis(rangeval = range(t_grid), nbasis = nbasis)
  Q_max <- numeric(R)
  for(b in 1:R){
    indices <- sample(1:n, n, replace = TRUE)
    X_boot <- simData@X[indices, ]
    fun_boot <- funData(argvals = list(t_grid), X = X_boot)
    fd_boot <- funData2fd(fun_boot, basisobj = basis_obj)
    pca_boot <- pca.fd(fd_boot, nharm = ncomp, centerfns = TRUE)
    scores <- pca_boot$scores
    coef_harm <- pca_boot$harmonics$coefs
    basis <- pca_boot$harmonics$basis
    mean_coef <- pca_boot$meanfd$coefs
    n_boot <- nrow(scores)
    K <- ncol(scores)
    gorros <- vector("list", n_boot)
    for(i in 1:n_boot){
      coef_i <- rep(0, nrow(coef_harm))
      for(k in 1:K){
        coef_i <- coef_i + scores[i,k] * coef_harm[,k]
      }
      coef_i <- coef_i + mean_coef
      gorros[[i]] <- fd(coef_i, basis)
    }
    coef_matrix <- do.call(cbind, lapply(gorros, function(f) f$coefs))
    gorros_fd <- fd(coef_matrix, basis)
    residuos <- fd_boot - gorros_fd
    Q_vals <- diag(inprod(residuos, residuos))
    Q_max[b] <- max(Q_vals)
  }
  quantile(Q_max, probs = 1 - alpha)
}

UCL_Q <- UCL_Q_FPCA_sim(simData = escenario1$simData, alpha = alpha, ncomp = 3, R = 3000)

# Carta Q
plot(Q, pch = 16, xlab = "Observación", ylab = "Q", main = "Carta Q", ylim=c(0,UCL_Q))
abline(h = UCL_Q, col = "red", lwd = 2)

Code 
# Proceso bajo control

# Índices de Capacidad

ncomp <- 3

# Base
basis_obj <- create.fourier.basis(c(0,1), nbasis = 20)

USL_fd_uni <- Data2fd(argvals = t, y = USL_uni, basisobj = basis_obj)
LSL_fd_uni <- Data2fd(argvals = t, y = LSL_uni, basisobj = basis_obj)
mean_fd <- Data2fd(argvals = t, y = mu, basisobj = basis_obj)

USL_centered_uni <- USL_fd_uni - mean_fd
LSL_centered_uni <- LSL_fd_uni - mean_fd

USL_pc_uni <- inprod(USL_centered_uni, pca$harmonics)
LSL_pc_uni <- inprod(LSL_centered_uni, pca$harmonics)

# Límites
lower_uni <- pmin(LSL_pc_uni, USL_pc_uni)
upper_uni <- pmax(LSL_pc_uni, USL_pc_uni)

# Índices de capacidad
cuantiles <- function(x, probs = c(0.005, 0.5, 0.995)){
  dens <- density(x, bw = "SJ")
  dx <- dens$x
  dy <- dens$y
  cdf <- cumsum(dy) * mean(diff(dx))
  cdf <- cdf / max(cdf)
  approx(cdf, dx, xout = probs)$y
}

# Cuantiles
cuantiles_proc <- matrix(NA, nrow = ncomp, ncol = 3)
for(i in 1:ncomp){
  cuantiles_proc[i,] <- cuantiles(scores[,i])
}

# Índices Cpk
indices_cpk_uni <- numeric(ncomp)
for(i in 1:ncomp){
  indices_cpk_uni[i] <- min((upper_uni[i] - cuantiles_proc[i,2]) / (cuantiles_proc[i,3] - cuantiles_proc[i,2]),
                            (cuantiles_proc[i,2] - lower_uni[i]) / (cuantiles_proc[i,2] - cuantiles_proc[i,1]))
}

MCnpk_uni <- prod(indices_cpk_uni)^(1/ncomp)
indices_cpk_uni
[1] 1.7475613 0.2328183 0.4310340
Code 
MCnpk_uni
[1] 0.5597409
Code 
# Gráfico especificaciones

matplot(t, t(escenario1$simData@X), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "")
lines(t, USL_uni, col = "red", lwd = 3)
lines(t, LSL_uni, col = "red", lwd = 3)
lines(t, colMeans(sim_uni$simData@X), lwd = 3, col = "blue")

Code 
# Gráfico especificaciones en espacio de scores

par(mfrow = c(1, ncomp))

for(m in 1:ncomp){
  plot(scores[,m], pch = 16, cex = 0.8, col = "black", xlab = "Observación",
       ylab = bquote(rho[.(m)]), main = paste("Score", m), ylim=c(min(min(scores), lower_uni[m]), max(max(scores), upper_uni[m]) ))
  abline(h = lower_uni[m], col = "blue", lwd = 2, lty = 2)
  abline(h = upper_uni[m], col = "red", lwd = 2, lty = 2)
}

Función para todos los escenarios

Code 
run_scenario_fpca <- function(M_sim = 3, eFunType = "Fourier", eValType = "exponential", basis_type = "Fourier",
                              N = 250, ncomp = 3, alpha = 0.01, R_boot = 3000, nbasis = 20, seed = 2143,
                              USL, LSL){
  
  set.seed(seed)
  t_grid <- seq(0,1,length.out = 100)
  sim <- simFunData(argvals = t_grid, M = M_sim, eFunType = eFunType, eValType = eValType, N = N)
  X <- sim$simData@X
  if(basis_type == "Fourier"){
    basis_obj <- create.fourier.basis(rangeval = c(0,1), nbasis = nbasis)
  } else{
    basis_obj <- create.bspline.basis(rangeval = c(0,1), nbasis = nbasis)
  }
  fd_obj <- funData2fd(sim$simData, basisobj = basis_obj)
  pca <- pca.fd(fd_obj, nharm = ncomp, centerfns = TRUE)
  scores <- pca$scores
  lambda <- pca$values
  T2 <- rowSums(sweep(scores[,1:ncomp]^2, 2, lambda[1:ncomp], "/"))
  alpha_t2 <- 1 - sqrt(1 - alpha)
  UCL_T2_FPCA_sim <- function(simData){
    n <- nObs(simData)
    T2_vmax <- numeric(R_boot)
    for(i in 1:R_boot){
      indices <- sample(1:n, n, replace = TRUE)
      X_boot <- simData@X[indices, ]
      fun_boot <- funData(argvals = t_grid, X = X_boot)
      fd_boot <- funData2fd(fun_boot, basisobj = basis_obj)
      pca_boot <- pca.fd( fd_boot, nharm = ncomp, centerfns = TRUE)
      scores_boot <- pca_boot$scores[,1:ncomp]
      lambda_boot <- pca_boot$values[1:ncomp]
      T2_vals <- rowSums(sweep(scores_boot^2, 2, lambda_boot, "/"))
      T2_vmax[i] <- max(T2_vals)
    }
    quantile(T2_vmax, probs = 1 - alpha_t2)
  }
  UCL_T2 <- UCL_T2_FPCA_sim(sim$simData)
  coef_harm <- pca$harmonics$coefs
  basis <- pca$harmonics$basis
  mean_coef <- pca$meanfd$coefs
  n <- nrow(scores)
  K <- ncol(scores)
  gorros <- vector("list", n)
  for(i in 1:n){
    coef_i <- rep(0, nrow(coef_harm))
    for(k in 1:K){
      coef_i <- coef_i + scores[i,k] * coef_harm[,k]
    }
    coef_i <- coef_i + mean_coef
    gorros[[i]] <- fd(coef_i, basis)
  }
  coef_matrix <- do.call(cbind, lapply(gorros, function(f) f$coefs))
  gorros_fd <- fd(coef_matrix, basis)
  residuos <- fd_obj - gorros_fd
  Q <- diag(inprod(residuos, residuos))
  UCL_Q_FPCA_sim <- function(simData){
    n <- nObs(simData)
    Q_max <- numeric(R_boot)
    for(b in 1:R_boot){
      indices <- sample(1:n, n, replace = TRUE)
      X_boot <- simData@X[indices, ]
      fun_boot <- funData(argvals = list(t_grid), X = X_boot)
      fd_boot <- funData2fd(fun_boot, basisobj = basis_obj)
      pca_boot <- pca.fd(fd_boot, nharm = ncomp, centerfns = TRUE)
      scores_boot <- pca_boot$scores
      coef_harm_boot <- pca_boot$harmonics$coefs
      basis_boot <- pca_boot$harmonics$basis
      mean_coef_boot <- pca_boot$meanfd$coefs
      n_boot <- nrow(scores_boot)
      K_boot <- ncol(scores_boot)
      gorros_boot <- vector("list", n_boot)
      for(i in 1:n_boot){
        coef_i <- rep(0, nrow(coef_harm_boot))
        for(k in 1:K_boot){
          coef_i <- coef_i + scores_boot[i,k] * coef_harm_boot[,k]
        }
        coef_i <- coef_i + mean_coef_boot
        gorros_boot[[i]] <- fd(coef_i, basis_boot)
      }
      coef_matrix_boot <- do.call(cbind, lapply(gorros_boot, function(f) f$coefs))
      gorros_fd_boot <- fd(coef_matrix_boot, basis_boot)
      residuos_boot <- fd_boot - gorros_fd_boot
      Q_vals <- diag(inprod(residuos_boot, residuos_boot))
      Q_max[b] <- max(Q_vals)
    }
    quantile(Q_max,
             probs = 1 - alpha)
  }
  UCL_Q <- UCL_Q_FPCA_sim(sim$simData)
  mu <- colMeans(X)
  USL_fd <- Data2fd(argvals = t_grid, y = USL, basisobj = basis_obj)
  LSL_fd <- Data2fd(argvals = t_grid, y = LSL, basisobj = basis_obj)
  mean_fd <- Data2fd(argvals = t_grid, y = mu, basisobj = basis_obj)
  USL_centered <- USL_fd - mean_fd
  LSL_centered <- LSL_fd - mean_fd
  USL_pc <- inprod(USL_centered, pca$harmonics)
  LSL_pc <- inprod(LSL_centered, pca$harmonics)
  lower <- pmin(LSL_pc, USL_pc)
  upper <- pmax(LSL_pc, USL_pc)
  cuantiles <- function(x, probs = c(0.005, 0.5, 0.995)){
    dens <- density(x, bw = "SJ")
    dx <- dens$x
    dy <- dens$y
    cdf <- cumsum(dy) * mean(diff(dx))
    cdf <- cdf / max(cdf)
    approx(cdf, dx, xout = probs)$y
  }
  cuantiles_proc <- matrix(NA, nrow = ncomp, ncol = 3)
  for(i in 1:ncomp){
    cuantiles_proc[i,] <- cuantiles(scores[,i])
  }
  
  indices_cpk <- numeric(ncomp)
  for(i in 1:ncomp){
    indices_cpk[i] <- min((upper[i] - cuantiles_proc[i,2]) / (cuantiles_proc[i,3] - cuantiles_proc[i,2]),
                          (cuantiles_proc[i,2] - lower[i]) / (cuantiles_proc[i,2] - cuantiles_proc[i,1]))
  }
  
  MCnpk <- prod(indices_cpk)^(1/ncomp)
  par(mfrow = c(1,1))
  plot(T2, pch = 16, xlab = "Observación", ylab = expression(T^2), main = expression("Carta " * T^2), ylim = c(0, max(T2, UCL_T2)))
  abline(h = UCL_T2, col = "red", lwd = 2,lty = 2)
  plot(Q, pch = 16, xlab = "Observación", ylab = "Q", main = "Carta Q", ylim=c(0, max(Q, UCL_Q)))
  abline(h = UCL_Q, col = "red", lwd = 2)
  
  # Curvas funcionales
  matplot(t_grid, t(X), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "X(t)", main = "Curvas funcionales")
  lines(t_grid, USL, col = "red", lwd = 3) 
  lines(t_grid, LSL, col = "blue", lwd = 3) 
  lines(t_grid, mu, col = "black", lwd = 3)
  
  # Scores
  par(mfrow = c(1,ncomp))
  for(m in 1:ncomp){
    plot(scores[,m], pch = 16, cex = 0.8, col = "black", xlab = "Observación",
         ylab = bquote(rho[.(m)]), main = paste("Score", m), ylim=c(min(min(scores), lower[m]), max(max(scores), upper[m])))
    abline(h = lower[m], col = "blue", lwd = 2, lty = 2)
    abline(h = upper[m], col = "red", lwd = 2, lty = 2)
  }
  
  return(list(sim = sim, pca = pca, T2 = T2, UCL_T2 = UCL_T2, Q = Q, 
              UCL_Q = UCL_Q, scores = scores, lower = lower, upper = upper,
              indices_cpk = indices_cpk, MCnpk = MCnpk))
}
Code 
# ESCENARIO 2: M = 5, eFunType = "Fourier", eValType = "exponential" --------

escenario2 <- run_scenario_fpca(M_sim = 5, eFunType = "Fourier", eValType = "exponential",
                                basis_type = "Fourier", ncomp = 3, USL = USL_uni, LSL = LSL_uni)

Code 
escenario2$indices_cpk
[1] 1.941541901 0.038055729 0.006979287
Code 
escenario2$MCnpk
[1] 0.08019106
Code 
# ESCENARIO 3: M = 5, eFunType = "Poly", eValType = "exponential" ---------

escenario3 <- run_scenario_fpca(M_sim = 5, eFunType = "Poly", eValType = "exponential", basis_type = "Fourier", ncomp = 3,
                                USL = USL_uni, LSL = LSL_uni)

Code 
escenario3$indices_cpk
[1] 1.94207107 0.03805423 0.18851709
Code 
escenario3$MCnpk
[1] 0.2406244
Code 
# ESCENARIO 4: M = 10, eFunType = "Fourier", eValType = "exponential" ----------------

escenario4 <- run_scenario_fpca(M_sim = 10, eFunType = "Fourier", eValType = "exponential", basis_type = "Fourier",
                                ncomp = 3, USL = USL_uni, LSL = LSL_uni)

Code 
escenario4$indices_cpk
[1] 1.76827877 0.58020738 0.08171738
Code 
escenario4$MCnpk
[1] 0.4376727
Code 
# ESCENARIO 5: M = 10, eFunType = "Poly", eValType = "exponential"

escenario5 <- run_scenario_fpca(M_sim = 10, eFunType = "Poly", eValType = "exponential", basis_type = "Fourier", ncomp = 3,
                                USL = USL_uni, LSL = LSL_uni)

Code 
escenario5$indices_cpk
[1] 1.7889992 0.4886865 0.2536972
Code 
escenario5$MCnpk
[1] 0.6053205
Code 
# ESCENARIO 6: M = 10, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline"

escenario6 <- run_scenario_fpca(M_sim = 10, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni, LSL = LSL_uni)

Code 
escenario6$indices_cpk
[1] 2.9280268 0.3347997 0.2707339
Code 
escenario6$MCnpk
[1] 0.6426398
Code 
# ESCENARIO 7: M = 10, eFunType = "Poly", eValType = "wiener", basis_type = "bspline" -------

escenario7 <- run_scenario_fpca(M_sim = 10, eFunType = "Poly", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni, LSL = LSL_uni)

Code 
escenario7$indices_cpk
[1] 2.9309130 0.1797672 0.3169782
Code 
escenario7$MCnpk
[1] 0.550699
Code 
# ESCENARIO 8: M = 30, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline" -------

escenario8 <- run_scenario_fpca(M_sim = 30, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni, LSL = LSL_uni)

Code 
escenario8$indices_cpk
[1] 3.28110564 0.09227339 0.38174027
Code 
escenario8$MCnpk
[1] 0.4871038
Code 
# ESCENARIO 9: M = 30, eFunType = "Poly", eValType = "wiener", basis_type = "bspline" -----------

escenario9 <- run_scenario_fpca(M_sim = 30, eFunType = "Poly", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni, LSL = LSL_uni)

Code 
escenario9$indices_cpk
[1]  3.28091050 -0.06221351  0.17669107
Code 
escenario9$MCnpk
[1] NaN
Code 
# ESCENARIO 10: M = 50, eFunType = "Fourier", eValType = "linear", basis_type = "bspline" --------------

escenario10 <- run_scenario_fpca(M_sim = 50, eFunType = "Fourier", eValType = "linear", basis_type = "bspline", ncomp = 3,
                                 nbasis = 9, USL = USL_uni, LSL = LSL_uni)

Code 
escenario10$indices_cpk
[1] 0.88162715 0.02376312 1.44096427
Code 
escenario10$MCnpk
[1] 0.3113727

Resultados

Code 
tabla_resultados_uni <- data.frame(
  
  Escenario = paste("Escenario", 1:10),
  
  M = c(3, 5, 5, 10, 10, 10, 10, 30, 30, 50),
  
  eFunType = c(
    "Fourier",
    "Fourier",
    "Poly",
    "Fourier",
    "Poly",
    "Fourier",
    "Poly",
    "Fourier",
    "Poly",
    "Fourier"
  ),
  
  eValType = c(
    "exponential",
    "exponential",
    "exponential",
    "exponential",
    "exponential",
    "wiener",
    "wiener",
    "wiener",
    "wiener",
    "linear"
  ),
  
  basis_type = c(
    "Fourier",
    "Fourier",
    "Fourier",
    "Fourier",
    "Fourier",
    "bspline",
    "bspline",
    "bspline",
    "bspline",
    "bspline"
  ),
  
  Cpk1 = c(
    indices_cpk_uni[1],
    escenario2$indices_cpk[1],
    escenario3$indices_cpk[1],
    escenario4$indices_cpk[1],
    escenario5$indices_cpk[1],
    escenario6$indices_cpk[1],
    escenario7$indices_cpk[1],
    escenario8$indices_cpk[1],
    escenario9$indices_cpk[1],
    escenario10$indices_cpk[1]
  ),
  
  Cpk2 = c(
    indices_cpk_uni[2],
    escenario2$indices_cpk[2],
    escenario3$indices_cpk[2],
    escenario4$indices_cpk[2],
    escenario5$indices_cpk[2],
    escenario6$indices_cpk[2],
    escenario7$indices_cpk[2],
    escenario8$indices_cpk[2],
    escenario9$indices_cpk[2],
    escenario10$indices_cpk[2]
  ),
  
  Cpk3 = c(
    indices_cpk_uni[3],
    escenario2$indices_cpk[3],
    escenario3$indices_cpk[3],
    escenario4$indices_cpk[3],
    escenario5$indices_cpk[3],
    escenario6$indices_cpk[3],
    escenario7$indices_cpk[3],
    escenario8$indices_cpk[3],
    escenario9$indices_cpk[3],
    escenario10$indices_cpk[3]
  ),
  
  MCnpk = c(
    MCnpk_uni,
    escenario2$MCnpk,
    escenario3$MCnpk,
    escenario4$MCnpk,
    escenario5$MCnpk,
    escenario6$MCnpk,
    escenario7$MCnpk,
    escenario8$MCnpk,
    escenario9$MCnpk,
    escenario10$MCnpk
  )
)

knitr::kable(tabla_resultados_uni)
Escenario M eFunType eValType basis_type Cpk1 Cpk2 Cpk3 MCnpk
Escenario 1 3 Fourier exponential Fourier 1.7475613 0.2328183 0.4310340 0.5597409
Escenario 2 5 Fourier exponential Fourier 1.9415419 0.0380557 0.0069793 0.0801911
Escenario 3 5 Poly exponential Fourier 1.9420711 0.0380542 0.1885171 0.2406244
Escenario 4 10 Fourier exponential Fourier 1.7682788 0.5802074 0.0817174 0.4376727
Escenario 5 10 Poly exponential Fourier 1.7889992 0.4886865 0.2536972 0.6053205
Escenario 6 10 Fourier wiener bspline 2.9280268 0.3347997 0.2707339 0.6426398
Escenario 7 10 Poly wiener bspline 2.9309130 0.1797672 0.3169782 0.5506990
Escenario 8 30 Fourier wiener bspline 3.2811056 0.0922734 0.3817403 0.4871038
Escenario 9 30 Poly wiener bspline 3.2809105 -0.0622135 0.1766911 NaN
Escenario 10 50 Fourier linear bspline 0.8816272 0.0237631 1.4409643 0.3113727

CASO MULTIVARIADO

En este caso se realizan varias simulaciones de observaciones funcionales para evaluar el desempeño y comportamiento del índice de capacidad multivariado basado en componentes principales dados los límites de especificación preestablecidos. En cada caso, se simulan \(N=250\) observaciones funcionales bivariadas tal que:

\[\boldsymbol{x}_i(\boldsymbol{t}) = \sum_{m=1}^M \rho_{im}\boldsymbol{\psi}_m(\boldsymbol{t})+\boldsymbol{\varepsilon}_i(\boldsymbol{t}), \quad \boldsymbol{\varepsilon}_i(\boldsymbol{t}) \overset{\text{iid}}{\sim} N_p(0,\sigma^2 \mathbf{I})\]

con \(\boldsymbol{t} \in T\) y \(i=1,...,250\).

Adicionalmente, los scores \(\rho_{im}\) son muestras independientes de una distribución \(N(0, \nu_m)\) para valores propios con decrecimiento lineal, exponencial o correspondientes a un proceso de Wiener.

Se consideran distintas configuraciones para cada escenario, donde varían:

  • Los valores de M
  • Tipo de funciones con que serán simuladas las curvas (“Poly”: Polinomios de Legendre o “Fourier”: Funciones de la base de Fourier).
  • El tipo de decrecimiento de los valores propios: (“linear”: \(\nu_m = \frac{M+1-m}{m}\), “exponential” = \(\nu_m = exp(-\frac{m-1}{2})\), “wiener” = \(\nu_m = \frac{1}{(\pi/2 \cdot (2m-1))^2}\)).
  • El tipo de base con que son ajustadas las observaciones simuladas (Fourier o B-Splines).

Límites de especificación

Se crea un conjunto de datos referencia para establecer los límites de especificación que serán utilizados en los escenarios propuestos de simulación. Para este caso, se consideran \(M=20\) funciones de la base de Fourier en el intervalo \([0,1]\) para cada uno de los procesos del proceso bivariado. Los límites de control están dados por:

\[USL(\boldsymbol{t})=\bar{\boldsymbol{x}}(\boldsymbol{t})+3.5 \boldsymbol{\sigma}(\boldsymbol{t})\] \[LSL(\boldsymbol{t})=\bar{\boldsymbol{x}}(\boldsymbol{t})-3.5 \boldsymbol{\sigma}(\boldsymbol{t})\]

Code 
# Construcción de límites de especificación

set.seed(2143)

t <- seq(0, 1, length.out = 100)

sim_ref <- simMultiFunData(type = "split", argvals = list(t,t), M = 20, eFunType = "Fourier",
                           eValType = "exponential", N = 250)

mu1 <- colMeans(sim_ref$simData[[1]]@X)
sd1 <- apply(sim_ref$simData[[1]]@X, 2, sd)

mu2 <- colMeans(sim_ref$simData[[2]]@X)
sd2 <- apply(sim_ref$simData[[2]]@X, 2, sd)

k <- 3.5

USL1 <- mu1 + k*sd1
LSL1 <- mu1 - k*sd1

USL2 <- mu2 + k*sd2
LSL2 <- mu2 - k*sd2

USL <- cbind(USL1, USL2)
LSL <- cbind(LSL1, LSL2)

par(mfrow=c(1,2))

for (j in 1:2){
  matplot(t, t(sim_ref$simData[[j]]@X), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "",
          ylim = c(-5, 5), main= paste("Proceso", j))
  
  lines(t, USL[,j], col = "red", lwd = 3)
  lines(t, LSL[,j], col = "red", lwd = 3)
  
  lines(t, colMeans(sim_ref$simData[[j]]@X), lwd = 3, col = "blue")
  
  legend("topright", legend = c("Límites de \n especificación"), col = c("red"), lwd = 3, lty = 1, bty = "n",
         cex = 0.4)
  
}

Code 
USL_fun <- multiFunData(list(funData(argvals = list(t), X = matrix(USL1, nrow = 1)),
                             funData(argvals = list(t), X = matrix(USL2, nrow = 1))))

LSL_fun <- multiFunData(list(funData(argvals = list(t), X = matrix(LSL1, nrow = 1)), 
                             funData(argvals = list(t), X = matrix(LSL2, nrow = 1))))

mean_list <- list(funData(argvals = list(t), X = matrix(mu1, nrow = 1)),
                  funData(argvals = list(t), X = matrix(mu2, nrow = 1)))

mean_mult <- multiFunData(mean_list)

En el primer escenario se considera \(M=3\), curvas simuladas con funciones de la base de Fourier, decrecimiento exponencial en los valores propios, curvas ajustadas con una base de Fourier.

Code 
# ESCENARIO 1: M = 3, type = "split, eValType = "exponential" ----------------------------------

set.seed(2143)

sim <- simMultiFunData(type = "split", argvals = list(seq(0,1,length.out=100),seq(0,1,length.out=100)), 
                       M = 3, eFunType = "Fourier", eValType = "exponential", N=250)

uni_exp <- vector("list", length = 2)
pca_uni <- vector("list", length = 2)
FourierBasis <- create.fourier.basis(rangeval = c(0,1), nbasis=20)

for (j in 1:2) {
  fd_j <- funData2fd(sim$simData[[j]], basisobj = FourierBasis)
  pca_j <- pca.fd(fd_j,nharm = 4,centerfns = TRUE)
  uni_exp[[j]] <- list(type = "given", 
                       functions = fd2funData(pca_j$harmonics, argvals = sim$simData[[j]]@argvals),
                       scores = pca_j$scores)
  pca_uni[[j]] <- pca_j
}

mfpca <- MFPCA(mFData = sim$simData, M=3, uniExpansions = uni_exp)
scores_mult <- mfpca$scores
lambda_mult <- mfpca$values
ncomp <- 3
T2_valores <- rowSums(sweep(scores_mult[,1:ncomp]^2, 2, lambda_mult[1:ncomp], "/"))

alpha = 1-sqrt(1-0.01)

# Límite de control
UCL_T2_MFPCA_sim <- function(simData, alpha = 0.01, ncomp = 5,
                             R = 2000, nbasis = 9){
  p <- length(simData)
  n <- nObs(simData)
  t_grid <- simData[[1]]@argvals[[1]]
  FourierBasis <- create.fourier.basis(
    rangeval = range(t_grid),
    nbasis = nbasis
  )
  
  T2_vmax <- numeric(R)
  
  for(i in 1:R){
    indices <- sample(1:n, n, replace = TRUE)
    lista_funData <- vector("list", p)
    for(k in 1:p){
      Xk <- simData[[k]]@X[indices, ]
      lista_funData[[k]] <- funData(argvals = list(t_grid), X = Xk)
    }
    mfData_boot <- multiFunData(lista_funData)
    uni_exp <- vector("list", p)
    for(j in 1:p){
      fd_j <- funData2fd(mfData_boot[[j]], basisobj = FourierBasis)
      pca_j <- pca.fd(fd_j, nharm = ncomp, centerfns = TRUE)
      uni_exp[[j]] <- list(
        type = "given",
        functions = fd2funData(pca_j$harmonics, argvals = mfData_boot[[j]]@argvals),
        scores = pca_j$scores
      )
    }
    mfpca <- MFPCA(
      mFData = mfData_boot,
      M = ncomp,
      uniExpansions = uni_exp
    )
    scores <- mfpca$scores[,1:ncomp]
    lambda <- mfpca$values[1:ncomp]
    T2_vals <- rowSums(sweep(scores^2, 2, lambda, "/"))
    T2_vmax[i] <- max(T2_vals)
  }
  quantile(T2_vmax, probs = 1 - alpha)
}

UCL_sim <- UCL_T2_MFPCA_sim(simData = sim$simData, alpha = alpha, 
                            ncomp = 3, R = 3000)

# Carta T2 

par(mfrow=c(1,1))
plot(T2_valores, type = "p", pch = 16, cex = 0.8, xlab = "Observación",
     ylab = expression(T^2), main = expression("Carta " * T^2),
     ylim=c(0,20))
abline(h = UCL_sim, col = "red", lwd = 2, lty = 2)

Code 
# Proceso bajo control


# Carta Q 

curvas_klmult <- predict(mfpca, scores = mfpca$scores)

p <- length(sim$simData)
n <- nObs(sim$simData)
t_grid <- sim$simData[[1]]@argvals[[1]]
dt <- diff(t_grid)[1]

Q_mult <- rep(0, n)
for(i in 1:p){
  X_orig <- sim$simData[[i]]@X
  X_hat <- curvas_klmult[[i]]@X
  error <- X_orig - X_hat
  Q_mult <- Q_mult + rowSums(error^2) * dt
}

UCL_Q_MFPCA_sim <- function(simData, alpha = 0.01, ncomp = 4,
                            R = 500, nbasis = 9){
  
  p <- length(simData)
  n <- nObs(simData)
  t_grid <- simData[[1]]@argvals[[1]]
  dt <- diff(t_grid)[1]
  
  FourierBasis <- create.fourier.basis(rangeval = range(t_grid),
                                       nbasis = nbasis)
  Q_max <- numeric(R)
  for(b in 1:R){
    indices <- sample(1:n, n, replace = TRUE)
    lista_funData <- vector("list", p)
    for(j in 1:p){
      Xj <- simData[[j]]@X[indices, ]
      lista_funData[[j]] <- funData(argvals = list(t_grid), X = Xj)
    }
    
    mfData_boot <- multiFunData(lista_funData)
    uni_exp <- vector("list", p)
    for(j in 1:p){
      fd_j <- funData2fd(mfData_boot[[j]], basisobj = FourierBasis)
      pca_j <- pca.fd(fd_j, nharm = ncomp, centerfns = TRUE)
      uni_exp[[j]] <- list(type = "given", functions = fd2funData(pca_j$harmonics, argvals = mfData_boot[[j]]@argvals),
                           scores = pca_j$scores)
    }

    mfpca_boot <- MFPCA(mFData = mfData_boot, M = ncomp, 
                        uniExpansions = uni_exp, fit = TRUE)
    
    residuos <- mfData_boot - mfpca_boot$fit
    Q_vals <- rep(0, n)
    for(j in 1:p){
      err_sq <- residuos[[j]]@X^2
      Q_vals <- Q_vals + rowSums(err_sq) * dt
    }
    Q_max[b] <- max(Q_vals)
  }
  quantile(Q_max, probs = 1 - alpha)
}

UCLQ_sim <- UCL_Q_MFPCA_sim(simData = sim$simData, alpha = alpha, 
                            ncomp = 3, R = 3000)

# Carta Q
plot(Q_mult, type = "p", pch = 16, cex = 0.6, xlab = "Observación",
     ylab = "Q", main = "Carta Q", ylim=c(0, UCLQ_sim+1.5))
abline(h = UCLQ_sim, col = "red", lwd = 2)

Code 
# Proceso bajo control

# Índices de Capacidad

ncomp <- 3

USL_centered <- USL_fun - mean_mult
LSL_centered <- LSL_fun - mean_mult

USL_pc <- scalarProduct(USL_centered, mfpca$functions)
LSL_pc <- scalarProduct(LSL_centered, mfpca$functions)

lower <- pmin(LSL_pc, USL_pc)
upper <- pmax(LSL_pc, USL_pc)

# Índice tipo Clements

cuantiles <- function(x, probs = c(0.005, 0.5, 0.995)) {
  dens <- density(x, bw = "SJ")
  dx <- dens$x
  dy <- dens$y
  cdf <- cumsum(dy) * mean(diff(dx))
  cdf <- cdf / max(cdf)
  approx(cdf, dx, xout = probs)$y
}

cuantiles_proc <- matrix(NA, nrow = ncomp, ncol = 3)
for(i in 1:ncomp){
  cuantiles_proc[i,] <- cuantiles(scores_mult[,i])
}

indices_cpk <- numeric(ncomp)

for(i in 1:ncomp){
  indices_cpk[i] <- min((upper[i] - cuantiles_proc[i,2]) / (cuantiles_proc[i,3] - cuantiles_proc[i,2]),
                        (cuantiles_proc[i,2] - lower[i]) / (cuantiles_proc[i,2] - cuantiles_proc[i,1]))
}

MCnpk <- prod(indices_cpk)^(1/ncomp)
indices_cpk
[1] 0.1897962 0.1318090 2.5452291
Code 
MCnpk
[1] 0.3993189
Code 
# Gráfico observaciones y especificaciones

par(mfrow = c(1,2))
for(j in 1:2){
  Xj <- sim$simData[[j]]@X
  matplot(t_grid, t(Xj), type = "l", lty = 1, col = "grey",
    xlab = "t", ylab = "", main = paste("Componente", j), ylim=c(-6,6))
  lines(t_grid, USL[,j], col = "red", lwd = 3)
  lines(t_grid, LSL[,j], col = "blue", lwd = 3)
  lines(t_grid, colMeans(Xj), lwd = 3, col = "black")
}

Code 
# Gráfico scores y especificaciones en el espacio de scores

par(mfrow = c(1, ncomp))

for(m in 1:ncomp){
  plot(scores_mult[,m], pch = 16, cex = 0.8, col = "black", xlab = "Observación",
    ylab = bquote(rho[.(m)]), main = paste("Score", m))
  abline(h = lower[m], col = "blue", lwd = 2, lty = 2)
  abline(h = upper[m], col = "red", lwd = 2, lty = 2)
}

Función para simular cada escenario

Code 
# Función para simular cada escenario ---------------------------------------------

run_scenario_mfpca <- function(M_sim = 2, type_sim = "split", eFunType = "Fourier", basis_type = "Fourier", 
                               evalType_sim = "exponential", N = 250, ncomp = 2, alpha = 0.01, 
                               R_boot = 3000, nbasis = 9, seed = 1234, USL = NULL, LSL = NULL){
  set.seed(seed)
  
  # Simulación
  t_grid <- seq(0,1,length.out = 100)
  sim <- simMultiFunData(type = type_sim, argvals = list(t_grid, t_grid), M = M_sim,
                         eFunType = eFunType, eValType = evalType_sim, N = N)
  
  # FPCA univariado
  if(basis_type == "Fourier"){
    basis_obj <- create.fourier.basis(rangeval = c(0,1), nbasis = nbasis)
  } 
  if(basis_type == "bspline"){
    basis_obj <- create.bspline.basis(rangeval = c(0,1), nbasis = nbasis)
  }
  uni_exp <- vector("list", length = 2)
  pca_uni <- vector("list", length = 2)
  
  for(j in 1:2){
    fd_j <- funData2fd(sim$simData[[j]], basisobj = basis_obj)
    pca_j <- pca.fd(fd_j, nharm = ncomp, centerfns = TRUE)
    uni_exp[[j]] <- list(type = "given", functions = fd2funData(pca_j$harmonics, argvals = sim$simData[[j]]@argvals),
                         scores = pca_j$scores)
    pca_uni[[j]] <- pca_j
  }
  
  # MFPCA
  mfpca <- MFPCA(mFData = sim$simData, M = ncomp, uniExpansions = uni_exp, fit = TRUE)
  scores_mult <- mfpca$scores
  lambda_mult <- mfpca$values
  
  # T2
  T2_valores <- rowSums(sweep(scores_mult[,1:ncomp]^2, 2, lambda_mult[1:ncomp], "/"))
  alpha_t2 <- 1 - sqrt(1 - alpha)
  UCL_T2_MFPCA_sim <- function(simData){
    p <- length(simData)
    n <- nObs(simData)
    T2_vmax <- numeric(R_boot)
    for(i in 1:R_boot){
      indices <- sample(1:n, n, replace = TRUE)
      lista_funData <- vector("list", p)
      for(k in 1:p){
        Xk <- simData[[k]]@X[indices, ]
        lista_funData[[k]] <- funData(argvals = list(t_grid), X = Xk)
      }
      mfData_boot <- multiFunData(lista_funData)
      uni_exp_boot <- vector("list", p)
      for(j in 1:p){
        fd_j <- funData2fd(mfData_boot[[j]], basisobj = basis_obj)
        pca_j <- pca.fd(fd_j, nharm = ncomp, centerfns = TRUE)
        uni_exp_boot[[j]] <- list(
          type = "given",
          functions = fd2funData(pca_j$harmonics, argvals = mfData_boot[[j]]@argvals),
          scores = pca_j$scores)
      }
      mfpca_boot <- MFPCA(mFData = mfData_boot, M = ncomp, uniExpansions = uni_exp_boot)
      scores_boot <- mfpca_boot$scores[,1:ncomp]
      lambda_boot <- mfpca_boot$values[1:ncomp]
      T2_vals <- rowSums(sweep(scores_boot^2, 2, lambda_boot, "/"))
      T2_vmax[i] <- max(T2_vals)
    }
    quantile(T2_vmax, probs = 1 - alpha_t2)
  }
  UCL_T2 <- UCL_T2_MFPCA_sim(sim$simData)
  
  # Carta Q
  p <- length(sim$simData)
  n <- nObs(sim$simData)
  dt <- diff(t_grid)[1]
  Q_mult <- rep(0, n)
  for(i in 1:p){
    X_orig <- sim$simData[[i]]@X
    X_hat <- mfpca$fit[[i]]@X
    error <- X_orig - X_hat
    Q_mult <- Q_mult + rowSums(error^2) * dt
  }
  
  UCL_Q_MFPCA_sim <- function(simData){
    Q_max <- numeric(R_boot)
    for(b in 1:R_boot){
      indices <- sample(1:n, n, replace = TRUE)
      lista_funData <- vector("list", p)
      for(j in 1:p){
        Xj <- simData[[j]]@X[indices, ]
        lista_funData[[j]] <- funData(argvals = list(t_grid), X = Xj)
      }
      mfData_boot <- multiFunData(lista_funData)
      uni_exp_boot <- vector("list", p)
      
      for(j in 1:p){
        fd_j <- funData2fd(mfData_boot[[j]], basisobj = basis_obj)
        pca_j <- pca.fd(fd_j, nharm = ncomp, centerfns = TRUE)
        uni_exp_boot[[j]] <- list(type = "given", functions = fd2funData(pca_j$harmonics, argvals = mfData_boot[[j]]@argvals),
                                  scores = pca_j$scores)
      }
      
      mfpca_boot <- MFPCA(mFData = mfData_boot, M = ncomp, uniExpansions = uni_exp_boot, fit = TRUE)
      residuos <- mfData_boot - mfpca_boot$fit
      Q_vals <- rep(0, n)
      for(j in 1:p){
        err_sq <- residuos[[j]]@X^2
        Q_vals <- Q_vals + rowSums(err_sq) * dt
      }
      Q_max[b] <- max(Q_vals)
    }
    quantile(Q_max, probs = 1 - alpha)
  }
  UCL_Q <- UCL_Q_MFPCA_sim(sim$simData)
  
  # Especificaciones funcionales
  X1 <- sim$simData[[1]]@X
  X2 <- sim$simData[[2]]@X
  
  mu1 <- colMeans(X1)
  mu2 <- colMeans(X2)
  
  if(is.null(USL) | is.null(LSL)){
    stop("Especificar USL y LSL")
  }
  
  # Media multivariada
  mean_mult <- multiFunData(list(funData(argvals = list(t_grid), X = matrix(mu1, nrow = 1)),
                                 funData(argvals = list(t_grid), X = matrix(mu2, nrow = 1))))
  
  # Límites en el espacio de scores
  USL1 <- USL[ ,1]
  USL2 <- USL[ ,2]
  LSL1 <- LSL[ ,1]
  LSL2 <- LSL[ ,2]
  
  USL_fun <- multiFunData(list(funData(argvals = list(t_grid), X = matrix(USL1, nrow = 1)),
                               funData(argvals = list(t_grid), X = matrix(USL2, nrow = 1))))
  LSL_fun <- multiFunData(list(funData(argvals = list(t_grid), X = matrix(LSL1, nrow = 1)),
                               funData(argvals = list(t_grid), X = matrix(LSL2, nrow = 1))))
  
  USL_pc <- scalarProduct(USL_fun - mean_mult, mfpca$functions)
  LSL_pc <- scalarProduct(LSL_fun - mean_mult, mfpca$functions)
  
  lower <- pmin(LSL_pc, USL_pc)
  upper <- pmax(LSL_pc, USL_pc)
  
  # Índice de capacidad
  cuantiles <- function(x,probs = c(0.005, 0.5, 0.995)){
    dens <- density(x, bw = "SJ")
    dx <- dens$x
    dy <- dens$y
    cdf <- cumsum(dy) * mean(diff(dx))
    cdf <- cdf / max(cdf)
    approx(cdf, dx, xout = probs)$y
  }
  cuantiles_proc <- matrix(NA, nrow = ncomp, ncol = 3)
  for(i in 1:ncomp){
    cuantiles_proc[i,] <- cuantiles(scores_mult[,i])
  }
  indices_cpk <- numeric(ncomp)
  for(i in 1:ncomp){
    indices_cpk[i] <- min((upper[i] - cuantiles_proc[i,2]) / (cuantiles_proc[i,3] - cuantiles_proc[i,2]),
                          (cuantiles_proc[i,2] - lower[i]) / (cuantiles_proc[i,2] - cuantiles_proc[i,1]))
  }
  MCnpk <- prod(indices_cpk)^(1/ncomp)
  
  # Gráficos
  par(mfrow = c(1,2))
  for(j in 1:2){
    Xj <- sim$simData[[j]]@X
    matplot(t_grid, t(Xj), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "", main = paste("Componente", j))
    lines(t_grid, USL[,j], col = "red", lwd = 3)
    lines(t_grid, LSL[,j], col = "blue", lwd = 3)
    lines(t_grid, colMeans(Xj), col = "black", lwd = 3)
  }
  par(mfrow = c(1,1))
  plot(T2_valores, pch = 16, xlab = "Observación", ylab = expression(T^2), main = expression("Carta " * T^2),
       ylim = c(0,UCL_T2+5))
  abline(h = UCL_T2, col = "red", lwd = 2, lty = 2)
  
  plot(Q_mult, pch = 16, xlab = "Observación", ylab = "Q", main = "Carta Q", ylim = c(0, UCL_Q+1.5))
  abline(h = UCL_Q, col = "red", lwd = 2)
  
  par(mfrow = c(1,ncomp))
  for(m in 1:ncomp){
    plot(scores_mult[,m], pch = 16, xlab = "Observación", ylab = bquote(rho[.(m)]), main = paste("Score", m),
         ylim = c(min(min(scores_mult[,m]), lower[m])-5, max(max(scores_mult[,m]), upper[m])+5))
    abline(h = lower[m], col = "blue", lwd = 2, lty = 2)
    abline(h = upper[m], col = "red", lwd = 2, lty = 2)
  }
  
  # Resultados
  return(list(sim = sim, mfpca = mfpca, T2 = T2_valores, UCL_T2 = UCL_T2, Q = Q_mult, UCL_Q = UCL_Q, 
              scores_mult = scores_mult, indices_cpk = indices_cpk, MCnpk = MCnpk, lower = lower, 
              upper = upper))
}

Se considera \(M=5\), curvas simuladas con funciones de la base de Fourier, decrecimiento exponencial en los valores propios, curvas ajustadas con una base de Fourier.

Code 
# ESCENARIO 2: M=5, eFUnType = "Fourier", eValType = "exponential"-------------------

esc2 <- run_scenario_mfpca(M_sim = 5, type_sim = "split", evalType_sim = "exponential", 
                                 ncomp = 3, USL = USL, LSL = LSL)

Code 
esc2$indices_cpk
[1] 1.65095452 0.09626945 0.15640469
Code 
esc2$MCnpk
[1] 0.2918487

Se considera \(M=5\), curvas simuladas con funciones de Legendre, decrecimiento exponencial en los valores propios, curvas ajustadas con una base de Fourier.

Code 
# ESCENARIO 3: M=5, eFunType = "Poly", eValType = "exponential"-------------

esc3 <- run_scenario_mfpca(M_sim = 5, type_sim = "split", evalType_sim = "exponential", 
                                 ncomp = 3, USL = USL, LSL = LSL)

Code 
esc3$indices_cpk
[1] 1.65095452 0.09626945 0.15640469
Code 
esc3$MCnpk
[1] 0.2918487

Se considera \(M=10\), curvas simuladas con funciones de la base de Fourier, decrecimiento exponencial en los valores propios, curvas ajustadas con una base de Fourier.

Code 
# ESCENARIO 4: M=10, eFunType = "Fourier", eValType = "exponential" -------------

esc4 <- run_scenario_mfpca(M_sim = 10, type_sim = "split", evalType_sim = "exponential", 
                                 ncomp = 3, USL = USL, LSL = LSL)

Code 
esc4$indices_cpk
[1] 1.82309452 0.01691276 0.13491395
Code 
esc4$MCnpk
[1] 0.1608275

Se considera \(M=10\), curvas simuladas con polinomios de Legendre, decrecimiento exponencial en los valores propios, curvas ajustadas con una base de Fourier.

Code 
# ESCENARIO 5: M=10, eFunType = "Poly", eValType = "exponential"-------------

esc5 <- run_scenario_mfpca(M_sim = 10, type_sim = "split", eFunType = "Poly", evalType_sim = "exponential", 
                                 ncomp = 3, USL = USL, LSL = LSL)

Code 
esc5$indices_cpk
[1] 1.8111677 0.0961321 0.3223691
Code 
esc5$MCnpk
[1] 0.3828778

Se considera \(M=10\), curvas simuladas con funciones de la base de Fourier, valores propios de un proceso de Wiener, curvas ajustadas con una base de B-Splines.

Code 
# ESCENARIO 6: M=10, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline" ------------------

esc6 <- run_scenario_mfpca(M_sim = 10, type = "split", eFunType = "Fourier", 
                       basis_type = "bspline", evalType_sim = "wiener", 
                       ncomp = 3, alpha = 0.01, nbasis = 9, USL = USL, LSL = LSL)

Code 
esc6$indices_cpk
[1]  2.8326441  0.2864970 -0.1487736
Code 
esc6$MCnpk
[1] NaN

Se considera \(M=10\), curvas simuladas con polinomios de Legendre, valores propios de un proceso de Wiener, curvas ajustadas con una base B-Splines.

Code 
# ESCENARIO 7: M=10, eFunType = "Poly", eValType = "wiener", basis_type = "bspline" ------------------

esc7 <- run_scenario_mfpca(M_sim = 10, type = "split", eFunType = "Poly", 
                                 basis_type = "bspline", evalType_sim = "wiener", 
                                 ncomp = 3, alpha = 0.01, nbasis = 9, USL = USL, LSL = LSL)

Code 
esc7$indices_cpk
[1]  2.8283346 -0.1214782  0.8681101
Code 
esc7$MCnpk
[1] NaN

Se considera \(M=30\), curvas simuladas con funciones de la base de Fourier, valores propios de un proceso de Wiener, curvas ajustadas con una base B-Splines.

Code 
# ESCENARIO 8: M=30, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline" ------------------

esc8 <- run_scenario_mfpca(M_sim = 30, type = "split", eFunType = "Fourier", 
                                 basis_type = "bspline", evalType_sim = "wiener", 
                                 ncomp = 3, alpha = 0.01, nbasis = 9, USL = USL, LSL = LSL)

Code 
esc8$indices_cpk
[1] 3.1767340 0.0384601 0.2193286
Code 
esc8$MCnpk
[1] 0.2992464

Se considera \(M=30\), curvas simuladas con polinomios de Lengendre, valores propios de un proceso de Wiener, curvas ajustadas con una base B-Splines.

Code 
# ESCENARIO 9: M=30, eFunType = "Poly", eValType = "wiener", basis_type = "bspline" ------------------

esc9 <- run_scenario_mfpca(M_sim = 10, type = "split", eFunType = "Poly", 
                                 basis_type = "bspline", evalType_sim = "wiener", 
                                 ncomp = 3, alpha = 0.01, nbasis = 9, USL = USL, LSL = LSL)

Code 
esc9$indices_cpk
[1]  2.8283346 -0.1214782  0.8681101
Code 
esc9$MCnpk
[1] NaN

Se considera \(M=50\), curvas simuladas con funciones de la base de Fourier, decrecimiento lineal en los valores propios, curvas ajustadas con una base B-Splines.

Code 
# ESCENARIO 10: M=50, eFunType = "Fourier", eValType = "linear", basis_type = "bspline" ------------------

esc10 <- run_scenario_mfpca(M_sim = 50, type = "split", eFunType = "Fourier", 
                                 basis_type = "bspline", evalType_sim = "linear", 
                                 ncomp = 3, alpha = 0.01, nbasis = 9, USL = USL, LSL = LSL)

Code 
esc10$indices_cpk
[1] 0.08863175 1.32614508 0.31606771
Code 
esc10$MCnpk
[1] 0.3336723

Se considera \(M=50\), curvas simuladas con polinomios de Legendre, decrecimiento exponencial en los valores propios, curvas ajustadas con una base B-Splines.

Code 
# ESCENARIO 11: M=50, eFunType = "Poly", eValType = "exponential", basis_type = "bspline" ------------------

esc11 <- run_scenario_mfpca(M_sim = 10, type = "split", eFunType = "Fourier", 
                                  basis_type = "bspline", evalType_sim = "linear", 
                                  ncomp = 3, alpha = 0.01, nbasis = 9, USL = USL, LSL = LSL)

Code 
esc11$indices_cpk
[1] 1.83005814 0.42498544 0.09522876
Code 
esc11$MCnpk
[1] 0.4199546

Resultados

Code 
resultados_1 <- data.frame(
  
  Escenario = paste("Escenario", 1:11),
  
  M = c(3,5,5,10,10,10,10,30,30,50,50),
  
  eFunType = c("Fourier","Fourier","Poly","Fourier","Poly",
               "Fourier","Poly","Fourier","Poly",
               "Fourier","Poly"),
  
  eValType = c("exponential","exponential","exponential",
               "exponential","exponential",
               "wiener","wiener",
               "wiener","wiener",
               "linear","linear"),
  
  basis_type = c("Fourier","Fourier","Fourier",
                 "Fourier","Fourier",
                 "bspline","bspline",
                 "bspline","bspline",
                 "bspline","bspline"),
  
  Cpk_1 = c(indices_cpk[1],
    esc2$indices_cpk[1],
    esc3$indices_cpk[1],
    esc4$indices_cpk[1],
    esc5$indices_cpk[1],
    esc6$indices_cpk[1],
    esc7$indices_cpk[1],
    esc8$indices_cpk[1],
    esc9$indices_cpk[1],
    esc10$indices_cpk[1],
    esc11$indices_cpk[1]
  ),
  
  Cpk_2 = c(indices_cpk[2],
    esc2$indices_cpk[2],
    esc3$indices_cpk[2],
    esc4$indices_cpk[2],
    esc5$indices_cpk[2],
    esc6$indices_cpk[2],
    esc7$indices_cpk[2],
    esc8$indices_cpk[2],
    esc9$indices_cpk[2],
    esc10$indices_cpk[2],
    esc11$indices_cpk[2]
  ),
  
  Cpk_3 = c(indices_cpk[3],
    esc2$indices_cpk[3],
    esc3$indices_cpk[3],
    esc4$indices_cpk[3],
    esc5$indices_cpk[3],
    esc6$indices_cpk[3],
    esc7$indices_cpk[3],
    esc8$indices_cpk[3],
    esc9$indices_cpk[3],
    esc10$indices_cpk[3],
    esc11$indices_cpk[3]
  ),
  
  MCnpk = c(MCnpk,
    esc2$MCnpk,
    esc3$MCnpk,
    esc4$MCnpk,
    esc5$MCnpk,
    esc6$MCnpk,
    esc7$MCnpk,
    esc8$MCnpk,
    esc9$MCnpk,
    esc10$MCnpk,
    esc11$MCnpk
  )
)

resultados_1[,6:9] <- round(resultados_1[,6:9], 3)

knitr::kable(resultados_1)
Escenario M eFunType eValType basis_type Cpk_1 Cpk_2 Cpk_3 MCnpk
Escenario 1 3 Fourier exponential Fourier 0.190 0.132 2.545 0.399
Escenario 2 5 Fourier exponential Fourier 1.651 0.096 0.156 0.292
Escenario 3 5 Poly exponential Fourier 1.651 0.096 0.156 0.292
Escenario 4 10 Fourier exponential Fourier 1.823 0.017 0.135 0.161
Escenario 5 10 Poly exponential Fourier 1.811 0.096 0.322 0.383
Escenario 6 10 Fourier wiener bspline 2.833 0.286 -0.149 NaN
Escenario 7 10 Poly wiener bspline 2.828 -0.121 0.868 NaN
Escenario 8 30 Fourier wiener bspline 3.177 0.038 0.219 0.299
Escenario 9 30 Poly wiener bspline 2.828 -0.121 0.868 NaN
Escenario 10 50 Fourier linear bspline 0.089 1.326 0.316 0.334
Escenario 11 50 Poly linear bspline 1.830 0.425 0.095 0.420