Simulaciones

Author

Ana María Valencia

Librerías

library(fda)
library(MFPCA)
library(funData)
library(depthTools)
library(roahd)
library(Matrix)
library(MASS)
library(corrplot)

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:

\[x_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}\)).

Límites de especificación (utilizando la media)

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))

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 = 1000, 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))
}

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=5\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de Fourier. Los valores propios se generan con un decrecimiento exponencial.

Code 
# ESCENARIO 1: 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

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=5\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un decrecimiento exponencial.

Code 
# ESCENARIO 2: 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

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=10\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un decrecimiento exponencial.

Code 
# ESCENARIO 3: 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

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=10\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

Code 
# ESCENARIO 4: 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

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=30\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

Code 
# ESCENARIO 5: 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

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=30\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=50\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un decrecimiento lineal.

Code 
# ESCENARIO 7: 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:7),
  
  M = c(5, 5, 10, 10, 30, 30, 50),
  
  eFunType = c(
    "Fourier",
    "Poly",
    "Poly",
    "Fourier",
    "Fourier",
    "Poly",
    "Fourier"
  ),
  
  eValType = c(
    "exponential",
    "exponential",
    "exponential",
    "wiener",
    "wiener",
    "wiener",
    "linear"
  ),
  
  Cpk1 = c(
    escenario2$indices_cpk[1],
    escenario3$indices_cpk[1],
    escenario5$indices_cpk[1],
    escenario6$indices_cpk[1],
    escenario8$indices_cpk[1],
    escenario9$indices_cpk[1],
    escenario10$indices_cpk[1]
  ),
  
  Cpk2 = c(
    escenario2$indices_cpk[2],
    escenario3$indices_cpk[2],
    escenario5$indices_cpk[2],
    escenario6$indices_cpk[2],
    escenario8$indices_cpk[2],
    escenario9$indices_cpk[2],
    escenario10$indices_cpk[2]
  ),
  
  Cpk3 = c(
    escenario2$indices_cpk[3],
    escenario3$indices_cpk[3],
    escenario5$indices_cpk[3],
    escenario6$indices_cpk[3],
    escenario8$indices_cpk[3],
    escenario9$indices_cpk[3],
    escenario10$indices_cpk[3]
  ),
  
  MCnpk = c(
    escenario2$MCnpk,
    escenario3$MCnpk,
    escenario5$MCnpk,
    escenario6$MCnpk,
    escenario8$MCnpk,
    escenario9$MCnpk,
    escenario10$MCnpk
  )
)

knitr::kable(tabla_resultados_uni)
Escenario M eFunType eValType Cpk1 Cpk2 Cpk3 MCnpk
Escenario 1 5 Fourier exponential 1.9415419 0.0380557 0.0069793 0.0801911
Escenario 2 5 Poly exponential 1.9420711 0.0380542 0.1885171 0.2406244
Escenario 3 10 Poly exponential 1.7889992 0.4886865 0.2536972 0.6053205
Escenario 4 10 Fourier wiener 2.9280268 0.3347997 0.2707339 0.6426398
Escenario 5 30 Fourier wiener 3.2811056 0.0922734 0.3817403 0.4871038
Escenario 6 30 Poly wiener 3.2809105 -0.0622135 0.1766911 NaN
Escenario 7 50 Fourier linear 0.8816272 0.0237631 1.4409643 0.3113727

Límites de especificación (utilizando la mediana)

En este caso, se considera el mismo conjunto de datos de referencia utilizado anteriormente. Sin embargo, los límites de especificación que serán utilizados en los escenarios propuestos de simulación se definen como:

\[USL(t)=\hat{med}_n(t)+3.5 \sigma(t)\] \[LSL(t)=\hat{med}_n(t)-3.5\sigma(t)\]

donde \(\hat{med}(t)\) hace referencia a la función mediana muestral, que se obtiene como:

\[\hat{med}_n = \arg \max_{x \in \{x_1,\ldots,x_n\}} MBD_n(x)\] donde \(MBD\) es la medida de profundidad Modified Band Depth.

Code 
X <- sim_uni$simData@X
mbd <- MBD(X)

Code 
# Mediana
mediana_indice <- mbd$ordering[1]

# Función mediana muestral
mediana <- X[mediana_indice, ]

k <- 3.5

# Límites de especificación
USL_uni_med <- mediana + k * sigma
LSL_uni_med <- mediana - 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_med, col = "red", lwd = 3)
lines(t, LSL_uni_med, col = "red", lwd = 3)
lines(t, mediana, lwd = 3, col = "blue")

Code 
USL_fun_uni_med <- funData(argvals = list(t), X = matrix(USL_uni_med, nrow = 1))
LSL_fun_uni_med <- funData(argvals = list(t), X = matrix(LSL_uni_med, nrow = 1))
mean_fun_uni <- funData(argvals = list(t), X = matrix(mu, nrow = 1))

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=5\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de Fourier. Los valores propios se generan con un decrecimiento exponencial.

Code 
# ESCENARIO 1: M = 5, eFunType = "Fourier", eValType = "exponential" --------

escenario2_med <- run_scenario_fpca(M_sim = 5, eFunType = "Fourier", eValType = "exponential",
                                basis_type = "Fourier", ncomp = 3, USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario2_med$indices_cpk
[1]  1.83324845 -0.05382534 -0.16759576
Code 
escenario2_med$MCnpk
[1] 0.2547751

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=5\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un decrecimiento exponencial.

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

escenario3_med <- run_scenario_fpca(M_sim = 5, eFunType = "Poly", eValType = "exponential", basis_type = "Fourier", ncomp = 3,
                                USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario3_med$indices_cpk
[1]  1.83703323 -0.05049960  0.02077519
Code 
escenario3_med$MCnpk
[1] NaN

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=10\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un decrecimiento exponencial.

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

escenario5_med <- run_scenario_fpca(M_sim = 10, eFunType = "Poly", eValType = "exponential", basis_type = "Fourier", ncomp = 3,
                                USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario5_med$indices_cpk
[1] 1.8079773 0.3493033 0.1338326
Code 
escenario5_med$MCnpk
[1] 0.4388532

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=10\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

escenario6_med <- run_scenario_fpca(M_sim = 10, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario6_med$indices_cpk
[1]  2.9070125  0.1328767 -0.3128785
Code 
escenario6_med$MCnpk
[1] NaN

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=30\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

escenario8_med <- run_scenario_fpca(M_sim = 30, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario8_med$indices_cpk
[1]  3.0948833 -0.1271266 -0.2287873
Code 
escenario8_med$MCnpk
[1] 0.4481645

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=30\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

escenario9_med <- run_scenario_fpca(M_sim = 30, eFunType = "Poly", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario9_med$indices_cpk
[1]  3.1007166 -0.3609460 -0.3340221
Code 
escenario9_med$MCnpk
[1] 0.720377

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=50\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un decrecimiento lineal.

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

escenario10_med <- run_scenario_fpca(M_sim = 50, eFunType = "Fourier", eValType = "linear", basis_type = "bspline", ncomp = 3,
                                 nbasis = 9, USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario10_med$indices_cpk
[1] 0.84764124 0.03309566 1.48982868
Code 
escenario10_med$MCnpk
[1] 0.3470349

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=50\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un decrecimiento lineal.

Code 
# ESCENARIO 8: M = 20, eFunType = "Fourier", eValType = "exponential", basis_type = "bspline" --------------

escenario11_med <- run_scenario_fpca(M_sim = 20, eFunType = "Fourier", eValType = "exponential", basis_type = "bspline", ncomp = 3,
                                 nbasis = 9, USL = USL_uni_med, LSL = LSL_uni_med)

Code 
escenario11_med$indices_cpk
[1]  1.91376855 -0.05947656  0.02069108
Code 
escenario11_med$MCnpk
[1] NaN

Resultados

Code 
tabla_resultados_uni_med <- data.frame(
  
  Escenario = paste("Escenario", 1:8),
  
  M = c(5, 5, 10, 10, 30, 30, 50, 20),
  
  eFunType = c(
    "Fourier",
    "Poly",
    "Poly",
    "Fourier",
    "Fourier",
    "Poly",
    "Fourier",
    "Fourier"
  ),
  
  eValType = c(
    "exponential",
    "exponential",
    "exponential",
    "wiener",
    "wiener",
    "wiener",
    "linear",
    "exponential"
  ),
  
  Cpk1 = c(
    escenario2_med$indices_cpk[1],
    escenario3_med$indices_cpk[1],
    escenario5_med$indices_cpk[1],
    escenario6_med$indices_cpk[1],
    escenario8_med$indices_cpk[1],
    escenario9_med$indices_cpk[1],
    escenario10_med$indices_cpk[1],
    escenario11_med$indices_cpk[1]
  ),
  
  Cpk2 = c(
    escenario2_med$indices_cpk[2],
    escenario3_med$indices_cpk[2],
    escenario5_med$indices_cpk[2],
    escenario6_med$indices_cpk[2],
    escenario8_med$indices_cpk[2],
    escenario9_med$indices_cpk[2],
    escenario10_med$indices_cpk[2],
    escenario11_med$indices_cpk[2]
  ),
  
  Cpk3 = c(
    escenario2_med$indices_cpk[3],
    escenario3_med$indices_cpk[3],
    escenario5_med$indices_cpk[3],
    escenario6_med$indices_cpk[3],
    escenario8_med$indices_cpk[3],
    escenario9_med$indices_cpk[3],
    escenario10_med$indices_cpk[3],
    escenario11_med$indices_cpk[3]
  ),
  
  MCnpk = c(
    escenario2_med$MCnpk,
    escenario3_med$MCnpk,
    escenario5_med$MCnpk,
    escenario6_med$MCnpk,
    escenario8_med$MCnpk,
    escenario9_med$MCnpk,
    escenario10_med$MCnpk,
    escenario11_med$MCnpk
    
  )
)

knitr::kable(tabla_resultados_uni_med)
Escenario M eFunType eValType Cpk1 Cpk2 Cpk3 MCnpk
Escenario 1 5 Fourier exponential 1.8332485 -0.0538253 -0.1675958 0.2547751
Escenario 2 5 Poly exponential 1.8370332 -0.0504996 0.0207752 NaN
Escenario 3 10 Poly exponential 1.8079773 0.3493033 0.1338326 0.4388532
Escenario 4 10 Fourier wiener 2.9070125 0.1328767 -0.3128785 NaN
Escenario 5 30 Fourier wiener 3.0948833 -0.1271266 -0.2287873 0.4481645
Escenario 6 30 Poly wiener 3.1007166 -0.3609460 -0.3340221 0.7203770
Escenario 7 50 Fourier linear 0.8476412 0.0330957 1.4898287 0.3470349
Escenario 8 20 Fourier exponential 1.9137686 -0.0594766 0.0206911 NaN

Límites de especificación (utilizando cuantiles)

Code 
# Con la relación p=|2q-1|
p <- abs(2*0.975 - 1)
n_esp <- dim(sim_uni$simData@X)[1]
k1 <- ceiling(p * n_esp)
indices_esp <- mbd$ordering[1:k1]
curvas <- sim_uni$simData@X[indices_esp, ]
LSLq <- apply(curvas, 2, min)
USLq <- apply(curvas, 2, max)

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

Code 
USLq_fun <- funData(argvals = list(t), X = matrix(USLq, nrow = 1))
LSLq_fun <- funData(argvals = list(t), X = matrix(LSLq, nrow = 1))

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=5\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de Fourier. Los valores propios se generan con un decrecimiento exponencial.

Code 
# ESCENARIO 1: M = 5, eFunType = "Fourier", eValType = "exponential" --------

escenario2q <- run_scenario_fpca(M_sim = 5, eFunType = "Fourier", eValType = "exponential",
                                basis_type = "Fourier", ncomp = 3, USL = USLq, LSL = LSLq)

Code 
escenario2q$indices_cpk
[1] 1.51828550 0.05457561 0.05531880
Code 
escenario2q$MCnpk
[1] 0.1661148

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=5\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un decrecimiento exponencial.

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

escenario3q <- run_scenario_fpca(M_sim = 5, eFunType = "Poly", eValType = "exponential", basis_type = "Fourier", ncomp = 3,
                                USL = USLq, LSL = LSLq)

Code 
escenario3q$indices_cpk
[1] 1.51545659 0.13242535 0.05763783
Code 
escenario3q$MCnpk
[1] 0.2261556

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=10\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un decrecimiento exponencial.

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

escenario5q <- run_scenario_fpca(M_sim = 10, eFunType = "Poly", eValType = "exponential", basis_type = "Fourier", ncomp = 3,
                                USL = USLq, LSL = LSLq)

Code 
escenario5q$indices_cpk
[1] 1.4003394 0.2853110 0.1031544
Code 
escenario5q$MCnpk
[1] 0.3454193

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=10\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

escenario6q <- run_scenario_fpca(M_sim = 10, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USLq, LSL = LSLq)

Code 
escenario6q$indices_cpk
[1] 2.2682000 0.2644760 0.7422495
Code 
escenario6q$MCnpk
[1] 0.7636116

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=30\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

escenario8q <- run_scenario_fpca(M_sim = 30, eFunType = "Fourier", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USLq, LSL = LSLq)

Code 
escenario8q$indices_cpk
[1] 2.57279146 0.04341968 0.68395456
Code 
escenario8q$MCnpk
[1] 0.4243324

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=30\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de polinomios de Legendre. Los valores propios se generan con un comportamiento como el de un proceso de Wiener.

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

escenario9q <- run_scenario_fpca(M_sim = 30, eFunType = "Poly", eValType = "wiener", basis_type = "bspline", ncomp = 3,
                                nbasis = 9, USL = USLq, LSL = LSLq)

Code 
escenario9q$indices_cpk
[1]  2.5738703  0.3214234 -0.1244616
Code 
escenario9q$MCnpk
[1] NaN

Se simulan observaciones funcionales univariadas utilizando una expansión de Karhunen-Loève con \(M=50\) componentes principales funcionales. Las funciones propias del proceso son generadas a partir de funciones de la base Fourier. Los valores propios se generan con un decrecimiento lineal.

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

escenario10q <- run_scenario_fpca(M_sim = 50, eFunType = "Fourier", eValType = "linear", basis_type = "bspline", ncomp = 3,
                                 nbasis = 9, USL = USLq, LSL = LSLq)

Code 
escenario10q$indices_cpk
[1]  0.701842810 -0.002412972  1.087470946
Code 
escenario10q$MCnpk
[1] NaN

Resultados

Code 
tabla_resultados_uni_q <- data.frame(
  
  Escenario = paste("Escenario", 1:7),
  
  M = c(5, 5, 10, 10, 30, 30, 50),
  
  eFunType = c(
    "Fourier",
    "Poly",
    "Poly",
    "Fourier",
    "Fourier",
    "Poly",
    "Fourier"
  ),
  
  eValType = c(
    "exponential",
    "exponential",
    "exponential",
    "wiener",
    "wiener",
    "wiener",
    "linear"
  ),
  
  Cpk1 = c(
    escenario2q$indices_cpk[1],
    escenario3q$indices_cpk[1],
    escenario5q$indices_cpk[1],
    escenario6q$indices_cpk[1],
    escenario8q$indices_cpk[1],
    escenario9q$indices_cpk[1],
    escenario10q$indices_cpk[1]
  ),
  
  Cpk2 = c(
    escenario2q$indices_cpk[2],
    escenario3q$indices_cpk[2],
    escenario5q$indices_cpk[2],
    escenario6q$indices_cpk[2],
    escenario8q$indices_cpk[2],
    escenario9q$indices_cpk[2],
    escenario10q$indices_cpk[2]
  ),
  
  Cpk3 = c(
    escenario2q$indices_cpk[3],
    escenario3q$indices_cpk[3],
    escenario5q$indices_cpk[3],
    escenario6q$indices_cpk[3],
    escenario8q$indices_cpk[3],
    escenario9q$indices_cpk[3],
    escenario10q$indices_cpk[3]
  ),
  
  MCnpk = c(
    escenario2q$MCnpk,
    escenario3q$MCnpk,
    escenario5q$MCnpk,
    escenario6q$MCnpk,
    escenario8q$MCnpk,
    escenario9q$MCnpk,
    escenario10q$MCnpk
  )
)

knitr::kable(tabla_resultados_uni_q)
Escenario M eFunType eValType Cpk1 Cpk2 Cpk3 MCnpk
Escenario 1 5 Fourier exponential 1.5182855 0.0545756 0.0553188 0.1661148
Escenario 2 5 Poly exponential 1.5154566 0.1324254 0.0576378 0.2261556
Escenario 3 10 Poly exponential 1.4003394 0.2853110 0.1031544 0.3454193
Escenario 4 10 Fourier wiener 2.2682000 0.2644760 0.7422495 0.7636116
Escenario 5 30 Fourier wiener 2.5727915 0.0434197 0.6839546 0.4243324
Escenario 6 30 Poly wiener 2.5738703 0.3214234 -0.1244616 NaN
Escenario 7 50 Fourier linear 0.7018428 -0.0024130 1.0874709 NaN

CASO MULTIVARIADO

Simulación 1

En este caso se realizaron distintas simulaciones de procesos funcionales multivariados con el objetivo de evaluar el comportamiento del índice de capacidad multivariado basado en componentes principales considerando distintas estructuras de dependencia entre las componentes funcionales. En cada escenario se generan \(n=250\) observaciones funcionales bivariadas mediante la descomposición de Karhunen-Loève, es decir, las observaciones funcionales son de la forma:

\[ \boldsymbol{X}_i(t) = \left( X_i^{(1)}(t), X_i^{(2)}(t) \right), \qquad i=1,\dots,N \]

donde cada componente funcional se construye como:

\[ X_i^{(j)}(t) = \sum_{k=1}^{K} \xi_{ik}^{(j)}\phi_k(t), \qquad j=1,2 \] con \(t \in [0,1]\), \(\phi_k(t)\) funciones base ortogonales y \(\xi_{ik}^{(j)}\) los scores asociados a cada componente principal funcional.

Las funciones \(\phi_k(t)\) son tomadas de una base de Fourier, utilizada como aproximación de las eigenfunciones del proceso funcional.

La estructura de dependencia entre las componentes funcionales es controlada directamente mediante la simulación conjunta de los scores. Se define el vector:

\[ \boldsymbol{\xi}_i = ( \xi_{i1}^{(1)}, \dots, \xi_{iK}^{(1)}, \xi_{i1}^{(2)}, \dots, \xi_{iK}^{(2)} ) ^\top \]

y se asume que:

\[ \boldsymbol{\xi}_i \sim N(\boldsymbol{0},\Sigma) \]

donde \(\Sigma\) es una matriz de covarianza definida positiva con la cual se controla:

  • La variabilidad de cada componente principal funcional mediante los elementos de la diagonal.
  • La dependencia entre procesos funcionales mediante los bloques de covarianza cruzada.

En todos los escenarios, las covarianzas entre componentes principales pertenecientes al mismo proceso funcional se consideran nulas. Por lo tanto, la matriz \(\Sigma\) puede expresarse mediante una estructura en bloques de la forma:

\[ \Sigma = \begin{pmatrix} D_1 & C \\ C^\top & D_2 \end{pmatrix} \]

donde \(D_1\) y \(D_2\) son matrices diagonales asociadas a las varianzas de los scores de cada proceso funcional y \(C\) representa la dependencia entre procesos.

Cada componente funcional fue simulada utilizando la expansión de Karhunen-Loève para procesos funcionales univariados. Posteriormente, para cada componente funcional se realizó FPCA univariado considerando 3 componentes principales. A partir de estas expansiones univariadas se aplicó MFPCA, considerando 3 componentes principales multivariados.

Finalmente, los scores obtenidos fueron utilizados para la construcción de las cartas de control \(T^2\) y \(Q\), y para el cálculo del índice de capacidad multivariado basado en componentes principales funcionales.

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 = 200)

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)

Funciones utilizadas

Code 
# Función para simular observaciones -----------------------

sim_mfpca_fourier <- function(n = 100, nbasis = 3, M = 6, t = seq(0,1,length.out = 200),
                              Sigma = NULL, seed = NULL){
  if(!is.null(seed)){
    set.seed(seed)
  }
  
  if(is.null(Sigma)){
    Sigma <- matrix(c(
      # proceso 1
      2.0, 0.0, 0.0,  0.4, 0.1, 0.1,
      0.0, 1.0, 0.0,  0.2, 0.3, 0.2,
      0.0, 0.0, 0.5,  0.1, 0.2, 0.4,
      
      # proceso 2
      0.4, 0.2, 0.1,  2.0, 0.0, 0.0,
      0.1, 0.3, 0.2,  0.0, 1.0, 0.0,
      0.1, 0.2, 0.4,  0.0, 0.0, 0.5
    ), nrow = 6, byrow = TRUE)
  }
  
  eigvals <- eigen(Sigma)$values
  if(any(eigvals <= 0)){
    stop("Sigma no es definida positiva")
  }
  
  basis_obj <- create.fourier.basis(rangeval = c(0,1), nbasis = nbasis)
  Phi <- eval.basis(t, basis_obj)
  scores <- MASS::mvrnorm(n = n, mu = rep(0, 2*nbasis), Sigma = Sigma)

  Xi1 <- scores[,1:nbasis]
  Xi2 <- scores[,(nbasis+1):(2*nbasis)]
  
  X1 <- Xi1 %*% t(Phi)
  X2 <- Xi2 %*% t(Phi)
  
  fd1 <- funData(argvals = list(t), X = X1)
  fd2 <- funData(argvals = list(t), X = X2)
  mfd <- multiFunData(fd1, fd2)
  
  uni_exp <- vector("list", 2)
  pca_uni <- vector("list", 2)
  
  for(j in 1:2){
    fd_j <- funData2fd(mfd[[j]], basisobj = basis_obj)
    pca_j <- pca.fd(fd_j, nharm = nbasis, centerfns = TRUE)
    uni_exp[[j]] <- list(
      type = "given",
      functions = fd2funData(pca_j$harmonics, argvals = mfd[[j]]@argvals),
      scores = pca_j$scores
    )
    pca_uni[[j]] <- pca_j
  }

  fit_mfpca <- MFPCA(mFData = mfd, M = M, uniExpansions = uni_exp)
  reconstruction <- predict(fit_mfpca, scores = fit_mfpca$scores)
  
  return(list(t = t, Phi = Phi, Sigma = Sigma, scores = scores, Xi1 = Xi1,
              Xi2 = Xi2, X1 = X1, X2 = X2, mfd = mfd, pca_uni = pca_uni,
              fit_mfpca = fit_mfpca, reconstruction = reconstruction))
}



sim1 <- sim_mfpca_fourier(n = 100, nbasis = 3, M = 6, t = seq(0,1,length.out = 200),
                              seed = 123)

### Función para cartas de control ------------------------------------------

analisis_mfpca <- function(sim_obj, alpha = 0.01, ncomp = 3, M = 3, R = 500, 
                           nbasis_monitor = 9
                           ){
  simData <- sim_obj$mfd
  p <- length(simData)
  
  uni_exp <- vector("list", length = p)
  pca_uni <- vector("list", length = p)
  
  FourierBasis <- create.fourier.basis(rangeval = c(0,1), nbasis = nbasis_monitor)
  
  for(j in 1:p){
    fd_j <- funData2fd(simData[[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 = simData[[j]]@argvals),
                         scores = pca_j$scores)
    pca_uni[[j]] <- pca_j
  }
  mfpca <- MFPCA(mFData = simData, M = M, uniExpansions = uni_exp)
  scores_mult <- mfpca$scores
  lambda_mult <- mfpca$values
  
  T2_valores <- rowSums(sweep(scores_mult[,1:M]^2, 2, lambda_mult[1:M], "/"))
  alpha_boot <- 1 - sqrt(1 - alpha)

  UCL_T2_MFPCA_sim <- function(simData, alpha, ncomp = 3, M = M, R = R, 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_boot <- MFPCA(mFData = mfData_boot, M = M, uniExpansions = uni_exp)
      scores_boot <- mfpca_boot$scores[,1:M]
      lambda_boot <- mfpca_boot$values[1:M]
      T2_vals <- rowSums(sweep(scores_boot^2, 2, lambda_boot, "/"))
      T2_vmax[i] <- max(T2_vals)
    }
    quantile(T2_vmax, probs = 1 - alpha)
  }
  
  UCL_T2 <- UCL_T2_MFPCA_sim(simData = simData, alpha = alpha_boot, ncomp = ncomp, M = M,
                             R = R, nbasis = nbasis_monitor)
  curvas_klmult <- predict(mfpca, scores = mfpca$scores)
  n <- nObs(simData)
  t_grid <- simData[[1]]@argvals[[1]]
  dt <- diff(t_grid)[1]
  Q_mult <- rep(0, n)
  for(i in 1:p){
    X_orig <- 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, ncomp, M, R, nbasis){
    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 = M, 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)
  }
  UCL_Q <- UCL_Q_MFPCA_sim(simData = simData, alpha = alpha_boot, ncomp = ncomp, M = M,
                           R = R, nbasis = nbasis_monitor)
  
  return(list(mfpca = mfpca, scores_mult = scores_mult, lambda_mult = lambda_mult, pca_uni = pca_uni,
              T2 = T2_valores, UCL_T2 = UCL_T2, Q = Q_mult, UCL_Q = UCL_Q,
              reconstruction = curvas_klmult))
}

### Función para calcular los índices de capacidad ---------------------------

indices_capacidad_mfpca <- function(sim_obj, analisis, USL, LSL, M = 3){
  
  t_grid <- sim_obj$t
  
  mean1 <- colMeans(sim_obj$X1)
  mean2 <- colMeans(sim_obj$X2)
  mean_fd1 <- funData( argvals = list(t_grid), X = matrix(mean1, nrow = 1))
  mean_fd2 <- funData(argvals = list(t_grid), X = matrix(mean2, nrow = 1))
  mean_mult <- multiFunData(mean_fd1, mean_fd2)
  
  USL_fd1 <- funData(argvals = list(t_grid), X = matrix(USL[,1], nrow = 1))
  USL_fd2 <- funData(argvals = list(t_grid), X = matrix(USL[,2], nrow = 1))
  LSL_fd1 <- funData(argvals = list(t_grid), X = matrix(LSL[,1], nrow = 1))
  LSL_fd2 <- funData(argvals = list(t_grid),  X = matrix(LSL[,2], nrow = 1))
  
  USL_fun <- multiFunData(USL_fd1, USL_fd2)
  LSL_fun <- multiFunData( LSL_fd1, LSL_fd2)
  
  USL_centered <- USL_fun - mean_mult
  LSL_centered <- LSL_fun - mean_mult
  
  mfpca <- analisis$mfpca
  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)
  
  scores_mult <- analisis$scores_mult
  
  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 = M, ncol = 3)
  for(i in 1:M){
    cuantiles_proc[i, ] <- cuantiles(scores_mult[,i])
  }
  
  indices_cpk <- numeric(M)
  for(i in 1:M){
    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/M)
  return(list(lower = lower, upper = upper, indices_cpk = indices_cpk, MCnpk = MCnpk))
}

### Función para graficar los scores ---------------------

graf_scores <- function(analisis, capacidad, M = 3){
  scores_mult <- analisis$scores_mult
  lower <- capacidad$lower
  upper <- capacidad$upper
  par(mfrow = c(1, M))
  for(m in 1:M){
    plot(scores_mult[,m], pch = 16, cex = 0.8,col = "black",
         xlab = "Observación", ylab = bquote(rho[.(m)]),
         main = paste("Score", m), ylim = c(min(min(scores_mult[,m]),lower[m]), max(max(scores_mult[,m]),upper[m])))
    abline(h = lower[m], col = "blue", lwd = 2, lty = 2)
    abline(h = upper[m], col = "red", lwd = 2, lty = 2)
  }
}

### Funcion para graficar los datos simulados con las especificaciones

graf_sim <- function(sim_obj, USL = NULL, LSL = NULL, ylim = NULL){
  t_grid <- sim_obj$t
  X_list <- list(sim_obj$X1, sim_obj$X2)
  p <- length(X_list)
  if(is.null(ylim)){
    ymin_curvas <- min(sapply(X_list, min))
    ymax_curvas <- max(sapply(X_list, max))
    if(!is.null(LSL)){
      ymin <- min(ymin_curvas, min(LSL))
    } else{
      ymin <- ymin_curvas
    }
    if(!is.null(USL)){
      ymax <- max(ymax_curvas, max(USL))
    } else{
      ymax <- ymax_curvas
    }
    ylim <- c(ymin, ymax)
  }
  par(mfrow = c(1,p))
  for(j in 1:p){
    Xj <- X_list[[j]]
    matplot(t_grid, t(Xj), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "",
            main = paste("Componente", j), ylim = ylim)
    lines(t_grid, colMeans(Xj), lwd = 3, col = "black")
    if(!is.null(USL)){lines(t_grid, USL[,j], col = "red", lwd = 3)
    }
    if(!is.null(LSL)){lines(t_grid, LSL[,j], col = "blue", lwd = 3)
    }
  }
}

scores_corr <- function(scores){
  corr_mat <- cor(scores)
  corrplot(corr_mat, method = "color", addCoef.col = "black",
           number.cex = 0.8, tl.col = "black", tl.srt = 0, diag = TRUE
  )
}

Se simulan observaciones que presenten baja variabilidad y baja correlación entre componentes funcionales.

Code 
Sigma1 <- matrix(c(
  
  0.5, 0.0, 0.0,  0.08, 0.03, 0.01,
  0.0, 0.3, 0.0,  0.03, 0.07, 0.02,
  0.0, 0.0, 0.1,  0.01, 0.02, 0.04,
  
  0.08, 0.03, 0.01,  0.5, 0.0, 0.0,
  0.03, 0.07, 0.02,  0.0, 0.3, 0.0,
  0.01, 0.02, 0.04,  0.0, 0.0, 0.1
  
), nrow = 6, byrow = TRUE)

sim1 <- sim_mfpca_fourier(n = 250 , nbasis = 3, M = 6, seed = 123, Sigma = Sigma1)
graf_sim(sim_obj = sim1, USL = USL, LSL = LSL)

Code 
analisis1 <- analisis_mfpca(sim_obj = sim1, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis1$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis1$UCL_T2))
abline(h = analisis1$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis1$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis1$UCL_Q))
abline(h = analisis1$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis1$reconstruction)

Code 
t_grid <- sim1$t
capacidad1 <- indices_capacidad_mfpca(sim_obj = sim1, analisis = analisis1,  
                                     USL = USL, LSL = LSL, M = 3)

capacidad1$indices_cpk
[1] 2.4213253 0.9354403 0.5676756
Code 
capacidad1$MCnpk
[1] 1.087401
Code 
graf_scores(analisis1, capacidad1, M=3)

Code 
scores_corr(analisis1$scores_mult)

# scores univariados
scores_uni1 <- cbind(analisis1$pca_uni[[1]]$scores, analisis1$pca_uni[[2]]$scores)
corrplot(cor(scores_uni1), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Se simulan observaciones que presenten variabilidad media y baja correlación entre componentes funcionales.

Code 
Sigma2 <- matrix(c(
  
  0.8, 0.0, 0.0,  0.25, 0.12, 0.05,
  0.0, 0.4, 0.0,  0.12, 0.20, 0.08,
  0.0, 0.0, 0.2,  0.05, 0.08, 0.12,
  
  0.25, 0.12, 0.05,  0.8, 0.0, 0.0,
  0.12, 0.20, 0.08,  0.0, 0.4, 0.0,
  0.05, 0.08, 0.12,  0.0, 0.0, 0.2
  
), nrow = 6, byrow = TRUE)

sim2 <- sim_mfpca_fourier(n = 250, nbasis = 3, M = 6, seed = 123, Sigma = Sigma2)
graf_sim(sim_obj = sim2, USL = USL, LSL = LSL)

Code 
analisis2 <- analisis_mfpca(sim_obj = sim2, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis2$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis2$UCL_T2))
abline(h = analisis2$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis2$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis2$UCL_Q))
abline(h = analisis2$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis2$reconstruction)

Code 
t_grid <- sim2$t
capacidad2 <- indices_capacidad_mfpca(sim_obj = sim2, analisis = analisis2,  
                                      USL = USL, LSL = LSL, M = 3)

capacidad2$indices_cpk
[1] 1.86165553 0.08132553 0.71820595
Code 
capacidad2$MCnpk
[1] 0.4773003
Code 
graf_scores(analisis2, capacidad2, M=3)

Code 
scores_corr(analisis2$scores_mult)

# scores univariados
scores_uni2 <- cbind(analisis2$pca_uni[[1]]$scores, analisis2$pca_uni[[2]]$scores)
corrplot(cor(scores_uni2), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Se simulan observaciones que presenten variabilidad moderada-alta y correlación moderada entre componentes funcionales.

Code 
Sigma3 <- matrix(c(
  
  1.0, 0.0, 0.0,   0.40, 0.22, 0.10,
  0.0, 0.6, 0.0,   0.22, 0.30, 0.12,
  0.0, 0.0, 0.3,   0.10, 0.12, 0.10,
  
  0.40, 0.22, 0.10,  1.2, 0.0, 0.0,
  0.22, 0.30, 0.12,  0.0, 0.64, 0.0,
  0.10, 0.12, 0.10,  0.0, 0.0, 0.2
  
), nrow = 6, byrow = TRUE)

sim3 <- sim_mfpca_fourier(n = 250, nbasis = 3, M = 6, seed = 123, Sigma = Sigma3)
graf_sim(sim_obj = sim3, USL = USL, LSL = LSL)

Code 
analisis3 <- analisis_mfpca(sim_obj = sim3, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis3$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis3$UCL_T2))
abline(h = analisis3$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis3$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis3$UCL_Q))
abline(h = analisis3$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis3$reconstruction)

Code 
t_grid <- sim3$t
capacidad3 <- indices_capacidad_mfpca(sim_obj = sim3, analisis = analisis3,  
                                      USL = USL, LSL = LSL, M = 3)

capacidad3$indices_cpk
[1] 1.4392968 0.5894751 0.8043269
Code 
capacidad3$MCnpk
[1] 0.8804056
Code 
graf_scores(analisis3, capacidad3, M=3)

Code 
scores_corr(analisis3$scores_mult)

# scores univariados
scores_uni3 <- cbind(analisis3$pca_uni[[1]]$scores, analisis3$pca_uni[[2]]$scores)
corrplot(cor(scores_uni3), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Se simulan observaciones que presenten alta variabilidad y correlación moderada-alta entre componentes funcionales.

Code 
Sigma4 <- matrix(c(
  
  2.5, 0.0, 0.0,  0.02, 0.35, 0.15,
  0.0, 1, 0.0,  0.25, 0.05, 0.20,
  0.0, 0.0, 0.4,  0.15, 0.20, 0.04,
  
  0.02, 0.25, 0.15,  1.5, 0.0, 0.0,
  0.35, 0.05, 0.20,  0.0, 0.8, 0.0,
  0.15, 0.20, 0.04,  0.0, 0.0, 0.2
  
), nrow = 6, byrow = TRUE)

sim4 <- sim_mfpca_fourier(n = 250, nbasis = 3, M = 6, seed = 123, Sigma = Sigma4)

analisis4 <- analisis_mfpca(sim_obj = sim4, alpha = 0.01, ncomp = 3, M = 3, R = 1000)
graf_sim(sim_obj = sim4, USL = USL, LSL = LSL)

Code 
# Carta T2
plot(analisis4$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis4$UCL_T2))
abline(h = analisis4$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis4$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis4$UCL_Q))
abline(h = analisis4$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis4$reconstruction)

Code 
t_grid <- sim4$t
capacidad4 <- indices_capacidad_mfpca(sim_obj = sim4, analisis = analisis4,  
                                      USL = USL, LSL = LSL, M = 3)

capacidad4$indices_cpk
[1] 0.6827604 1.3134108 0.4021197
Code 
capacidad4$MCnpk
[1] 0.7117729
Code 
graf_scores(analisis4, capacidad4, M=3)

Code 
scores_corr(analisis4$scores_mult)

# scores univariados
scores_uni4 <- cbind(analisis4$pca_uni[[1]]$scores, analisis4$pca_uni[[2]]$scores)
corrplot(cor(scores_uni4), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Se simulan observaciones que presenten alta variabilidad y alta correlación entre componentes funcionales.

Code 
Sigma5 <- matrix(c(
  
  2.4, 0.0, 0.0,   0.70, 0.35, 0.15,
  0.0, 0.8, 0.0,   0.35, 0.45, 0.18,
  0.0, 0.0, 0.5,   0.15, 0.18, 0.18,
  
  0.70, 0.35, 0.15,   2.6, 0.0, 0.0,
  0.35, 0.45, 0.18,   0.0, 1.0, 0.0,
  0.15, 0.18, 0.18,   0.0, 0.0, 0.4
  
), nrow = 6, byrow = TRUE)

sim5 <- sim_mfpca_fourier(n = 250, nbasis = 3, M = 6, seed = 123, Sigma = Sigma5)
graf_sim(sim_obj = sim5, USL = USL, LSL = LSL)

Code 
analisis5 <- analisis_mfpca(sim_obj = sim5, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis5$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis5$UCL_T2))
abline(h = analisis5$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis5$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis5$UCL_Q))
abline(h = analisis5$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis5$reconstruction)

Code 
t_grid <- sim5$t
capacidad5 <- indices_capacidad_mfpca(sim_obj = sim5, analisis = analisis5,  
                                      USL = USL, LSL = LSL, M = 3)

capacidad5$indices_cpk
[1] 1.0596171 0.3468498 0.2982798
Code 
capacidad5$MCnpk
[1] 0.4785986
Code 
graf_scores(analisis5, capacidad5, M=3)

Code 
scores_corr(analisis5$scores_mult)

# scores univariados
scores_uni5 <- cbind(analisis5$pca_uni[[1]]$scores, analisis5$pca_uni[[2]]$scores)
corrplot(cor(scores_uni5), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Se simulan observaciones que presenten alta variabilidad y baja correlación entre componentes funcionales.

Code 
Sigma6 <- matrix(c(
  
  1.8, 0.0, 0.0,   0.18, 0.08, 0.04,
  0.0, 1.1, 0.0,   0.08, 0.15, 0.06,
  0.0, 0.0, 0.7,   0.04, 0.06, 0.10,
  
  0.18, 0.08, 0.04,   1.7, 0.0, 0.0,
  0.08, 0.15, 0.06,   0.0, 1.0, 0.0,
  0.04, 0.06, 0.10,   0.0, 0.0, 0.6
  
), nrow = 6, byrow = TRUE)

sim6 <- sim_mfpca_fourier(n = 250, nbasis = 3, M = 6, seed = 123, Sigma = Sigma6)
graf_sim(sim_obj = sim6, USL = USL, LSL = LSL)

Code 
analisis6 <- analisis_mfpca(sim_obj = sim6, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis6$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis6$UCL_T2))
abline(h = analisis6$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis6$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis6$UCL_Q))
abline(h = analisis6$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis6$reconstruction)

Code 
t_grid <- sim6$t
capacidad6 <- indices_capacidad_mfpca(sim_obj = sim6, analisis = analisis6,  
                                      USL = USL, LSL = LSL, M = 3)

capacidad6$indices_cpk
[1] 1.1131298 1.0488950 0.2497074
Code 
capacidad6$MCnpk
[1] 0.6630858
Code 
graf_scores(analisis6, capacidad6, M=3)

Code 
scores_corr(analisis6$scores_mult)

# scores univariados
scores_uni6 <- cbind(analisis6$pca_uni[[1]]$scores, analisis6$pca_uni[[2]]$scores)
corrplot(cor(scores_uni6), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Resultados

Code 
resultados_dependencia <- data.frame(
  
  Escenario = c(
    "Escenario 1",
    "Escenario 2",
    "Escenario 3",
    "Escenario 4",
    "Escenario 5",
    "Escenario 6"
  ),
  
  Variabilidad = c(
    "Baja",
    "Media",
    "Moderada-Alta",
    "Alta",
    "Alta",
    "Alta"
  ),
  
  Covarianza = c(
    "Baja",
    "Moderada",
    "Moderada",
    "Moderada-Alta",
    "Alta",
    "Baja"
  ),
  
  Cpk_1 = c(
    capacidad1$indices_cpk[1],
    capacidad2$indices_cpk[1],
    capacidad3$indices_cpk[1],
    capacidad4$indices_cpk[1],
    capacidad5$indices_cpk[1],
    capacidad6$indices_cpk[1]
    
  ),
  
  Cpk_2 = c(
    capacidad1$indices_cpk[2],
    capacidad2$indices_cpk[2],
    capacidad3$indices_cpk[2],
    capacidad4$indices_cpk[2],
    capacidad5$indices_cpk[2],
    capacidad6$indices_cpk[2]
  ),
  
  Cpk_3 = c(
    capacidad1$indices_cpk[3],
    capacidad2$indices_cpk[3],
    capacidad3$indices_cpk[3],
    capacidad4$indices_cpk[3],
    capacidad5$indices_cpk[3],
    capacidad6$indices_cpk[3]
  ),
  
  MCpk = c(
    capacidad1$MCnpk,
    capacidad2$MCnpk,
    capacidad3$MCnpk,
    capacidad4$MCnpk,
    capacidad5$MCnpk,
    capacidad6$MCnpk
  )
)

resultados_dependencia[,4:7] <- round(
  resultados_dependencia[,4:7],
  3
)

knitr::kable(
  resultados_dependencia,
  caption = "Índices de capacidad para distintos escenarios de dependencia"
)
Índices de capacidad para distintos escenarios de dependencia
Escenario Variabilidad Covarianza Cpk_1 Cpk_2 Cpk_3 MCpk
Escenario 1 Baja Baja 2.421 0.935 0.568 1.087
Escenario 2 Media Moderada 1.862 0.081 0.718 0.477
Escenario 3 Moderada-Alta Moderada 1.439 0.589 0.804 0.880
Escenario 4 Alta Moderada-Alta 0.683 1.313 0.402 0.712
Escenario 5 Alta Alta 1.060 0.347 0.298 0.479
Escenario 6 Alta Baja 1.113 1.049 0.250 0.663

Simulación 2

Para simular las observaciones funcionales multivariadas, en primer lugar se generan las bases de funciones propias univariadas para cada uno de los \(p\) procesos funcionales. Para el \(j\)-ésimo proceso, se consideran \(M_j\) funciones base \(\phi_m^{(j)}\), las cuales representan las principales direcciones de variación de cada componente funcional.

Posteriormente, se define una matriz de correlaciones entre procesos

\[ R = \begin{pmatrix} 1 & r_{12} & \cdots & r_{1p} \\ r_{21} & 1 & \cdots & r_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ r_{p1} & r_{p2} & \cdots & 1 \end{pmatrix}, \]

donde cada elemento \(r_{ij}\) representa la correlación entre los procesos \(i\) y \(j\). Esta matriz se utiliza para construir una matriz de covarianzas conjunta \(\mathbf{Z}\) entre los scores univariados de todos los procesos.

Las entradas de \(\mathbf{Z}\) están dadas por

\[ Z_{mn}^{(ij)} = \operatorname{Cov}(\xi_m^{(i)},\xi_n^{(j)}). \]

Dentro de un mismo proceso, las covarianzas corresponden a los valores propios asociados a las funciones base univariadas. Para procesos distintos, las covarianzas se definen mediante

\[ Z_{mn}^{(ij)} = \text{Cov}(\xi_{m}^{(i)}, \xi_{n}^{(j)}) = r_{ij}\exp(-\alpha |m-n|) \sqrt{\lambda_m^{(i)}\lambda_n^{(j)}}, \]

donde \(\alpha > 0\) controla el decrecimiento de la correlación entre componentes de diferente orden. De esta forma, la correlación entre procesos se incorpora en la estructura de covarianza de los scores.

A partir de la matriz \(\mathbf{Z}\), se generan muestras de scores univariados utilizando una distribución normal multivariada. Luego, se realiza la descomposición espectral de \(\mathbf{Z}\) para obtener sus valores propios y vectores propios, los cuales permiten construir tanto los scores multivariados \(\rho_{im}\) como las funciones propias multivariadas \(\boldsymbol{\psi}_m\).

Finalmente, las observaciones funcionales multivariadas se generan utilizando la expansión truncada de Karhunen–Loève

\[ \mathbf{X}_i(\mathbf{t}) = \sum_{m=1}^{M} \rho_{im}\boldsymbol{\psi}_m(\mathbf{t}), \qquad i=1,\ldots,N. \]

Equivalentemente, para cada proceso \(j=1,\ldots,p\),

\[ X_i^{(j)}(t_j) = \sum_{m=1}^{M} \rho_{im}\psi_m^{(j)}(t_j). \]

Funciones utilizadas

Code 
MultiScores <- function(scores, Z, M){
  scores <- scale(scores, center = TRUE, scale = FALSE)
  eig <- eigen(Z, symmetric = TRUE)
  TrueVals <- eig$values[1:M]
  Cvec <- eig$vectors[, 1:M, drop = FALSE]
  rho <- scores %*% Cvec
  return(list(MultiScores = rho,
              TrueVals = TrueVals,
              Cvec = Cvec))
  
}

simMultiScores <- function(N, M, multiM, correlations, type){
  p <- length(M)
  Mtotal <- sum(M)
  # Matriz de correlaciones entre los procesos
  if (length(correlations) == 1) {
    R <- matrix(correlations, p, p)
    diag(R) <- 1
  } else {
    R <- matrix(1, p, p)
    R[upper.tri(R)] <- correlations
    R[lower.tri(R)] <- t(R)[lower.tri(R)]
  }
  # Valores propios
  lambdas <- lapply(M, function(m) {
    eVal(M = m, type = type)
  })
  # Construcción matriz de covarianzas
  bloques <- split(1:Mtotal, rep(1:p, M))
  Z <- matrix(0, Mtotal, Mtotal)
  ## Dentro de cada proceso
  for (j in 1:p) {
    Z[bloques[[j]], bloques[[j]]] <- diag(lambdas[[j]])
  }
  ## Covarianzas entre procesos
  for (j in 1:(p - 1)) {
    for (k in (j + 1):p) {
      for (m in 1:M[j]) {
        for (n in 1:M[k]) {
          a <- bloques[[j]][m]
          b <- bloques[[k]][n]
          rho_mn <- R[j, k] * exp(- 3 * abs(m - n))
          Z[a, b] <- rho_mn * sqrt(lambdas[[j]][m] * lambdas[[k]][n])
          Z[b, a] <- Z[a, b]
        }
      }
    }
  }
  if (min(eigen(Z, symmetric = TRUE)$values) < -1e-8) {
    stop("La matriz Z no es semidefinida positiva")
  }
  # Simular matriz de scores Xi
  Xi <- mvrnorm(n = N, mu = rep(0, Mtotal), Sigma = Z)
  # Construir los scores multivariados
  multiScores <- MultiScores(scores = Xi, Z = Z, M = multiM)
  return(list(UniScores = Xi,
              MultiScores = multiScores$MultiScores,
              TrueVals = multiScores$TrueVals,
              Cvec = multiScores$Cvec,
              Z = Z))
  
}

simUniBasis <- function(argvals, Ms, eFunType, eValType){
  p <- length(argvals)
  uniBasis <- vector("list", p)
  for(j in seq_len(p)){
    M <- if(is.list(Ms)) Ms[[j]] else Ms[j]
    eFun <- if(is.list(eFunType)) eFunType[[j]] else eFunType
    uniBasis[[j]] <- simFunData(N = 1, 
                                argvals = argvals[[j]], 
                                M = M, 
                                eFunType = eFun,
                                eValType = eValType)$trueFuns
  }
  return(uniBasis)
}

buildMultiFuns <- function(uniBasis, Cvec, M){
  p <- length(uniBasis)
  # Número de bases univariadas por proceso
  Ms <- sapply(uniBasis, nObs)
  # Bloques del vector propio asociados a cada proceso
  bloques <- split(seq_len(sum(Ms)), rep(seq_len(p), Ms))
  trueFuns <- vector("list", p)
  for(j in seq_len(p)){
    Phi <- uniBasis[[j]]@X
    dims_phi <- dim(Phi)
    # Matriz de funciones base
    Phi_mat <- matrix(Phi, nrow = Ms[j])
    # Coeficientes del proceso j en los vectores propios de Z
    C <- Cvec[bloques[[j]], 1:M, drop = FALSE]
    # Construcción de la función propia
    Psi_mat <- t(C) %*% Phi_mat
    # Objeto funData
    Psi <- array(Psi_mat, dim = c(M, dims_phi[-1]))
    trueFuns[[j]] <- funData(
      argvals = uniBasis[[j]]@argvals,
      X = Psi
    )
  }
  return(multiFunData(trueFuns))
}

simMultiCorr <- function(type, 
                         argvals, 
                         M,
                         Ms,
                         eFunType, 
                         eValType,
                         N,
                         correlations,
                         ignoreDeg = NULL){
  # Bases de funciones univariadas por proceso
  uniBasis <- simUniBasis(argvals = argvals, 
                          Ms = Ms, 
                          eFunType = eFunType, 
                          eValType = eValType)
  Mscores <- sapply(uniBasis, nObs)
  # Se generan los valores propios y scores con correlación 
  multiScores <- simMultiScores(N = N, 
                                M = Mscores, 
                                multiM = M,
                                correlations = correlations, 
                                type = eValType)
  trueVals <- multiScores$TrueVals
  scores <- multiScores$MultiScores
  Cvec <- multiScores$Cvec
  # Construcción de funciones propias multivariadas
  trueFuns <- buildMultiFuns(uniBasis = uniBasis, 
                             Cvec = Cvec,
                             M = M)
  Mtotal <- nObs(trueFuns)
  p <- length(trueFuns)
  # Simulación de los datos
  simData <- vector("list", p)
  for (j in seq_len(p)) {
    X <- apply(trueFuns[[j]]@X, -1, function(v) {
      scores %*% v
    })
    if (N == 1) {
      dim(X) <- c(1, nObsPoints(trueFuns[[j]]))
    }
    simData[[j]] <- funData(trueFuns[[j]]@argvals, X)
  }
  return(list(simData = multiFunData(simData), trueFuns = trueFuns, 
              trueVals = trueVals))
}


### Función para cartas de control ------------------------------------------

analisis_mfpca2 <- function(sim_obj, alpha = 0.01, ncomp = 3, M = 3, R = 500, 
                           nbasis_monitor = 9
                           ){
  simData <- sim_obj$simData
  p <- length(simData)
  
  uni_exp <- vector("list", length = p)
  pca_uni <- vector("list", length = p)
  
  FourierBasis <- create.fourier.basis(rangeval = c(0,1), nbasis = nbasis_monitor)
  
  for(j in 1:p){
    fd_j <- funData2fd(simData[[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 = simData[[j]]@argvals),
                         scores = pca_j$scores)
    pca_uni[[j]] <- pca_j
  }
  mfpca <- MFPCA(mFData = simData, M = M, uniExpansions = uni_exp)
  scores_mult <- mfpca$scores
  lambda_mult <- mfpca$values
  
  T2_valores <- rowSums(sweep(scores_mult[,1:M]^2, 2, lambda_mult[1:M], "/"))
  alpha_boot <- 1 - sqrt(1 - alpha)

  UCL_T2_MFPCA_sim <- function(simData, alpha, ncomp = 3, M = M, R = R, 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_boot <- MFPCA(mFData = mfData_boot, M = M, uniExpansions = uni_exp)
      scores_boot <- mfpca_boot$scores[,1:M]
      lambda_boot <- mfpca_boot$values[1:M]
      T2_vals <- rowSums(sweep(scores_boot^2, 2, lambda_boot, "/"))
      T2_vmax[i] <- max(T2_vals)
    }
    quantile(T2_vmax, probs = 1 - alpha)
  }
  
  UCL_T2 <- UCL_T2_MFPCA_sim(simData = simData, alpha = alpha_boot, ncomp = ncomp, M = M,
                             R = R, nbasis = nbasis_monitor)
  curvas_klmult <- predict(mfpca, scores = mfpca$scores)
  n <- nObs(simData)
  t_grid <- simData[[1]]@argvals[[1]]
  dt <- diff(t_grid)[1]
  Q_mult <- rep(0, n)
  for(i in 1:p){
    X_orig <- 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, ncomp, M, R, nbasis){
    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 = M, 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)
  }
  UCL_Q <- UCL_Q_MFPCA_sim(simData = simData, alpha = alpha_boot, ncomp = ncomp, M = M,
                           R = R, nbasis = nbasis_monitor)
  
  return(list(mfpca = mfpca, scores_mult = scores_mult, lambda_mult = lambda_mult, pca_uni = pca_uni,
              T2 = T2_valores, UCL_T2 = UCL_T2, Q = Q_mult, UCL_Q = UCL_Q,
              reconstruction = curvas_klmult))
}

### Función para calcular los índices de capacidad ---------------------------

indices_capacidad_mfpca2 <- function(sim_obj, analisis, USL, LSL, M = 3){

  simData <- sim_obj$simData
  p <- length(simData)

  mean_list <- vector("list", p)
  for(j in seq_len(p)){
    mean_j <- colMeans(simData[[j]]@X)
    mean_list[[j]] <- funData(
      argvals = simData[[j]]@argvals,
      X = matrix(mean_j, nrow = 1)
    )
  }
  mean_mult <- multiFunData(mean_list)
  USL_list <- vector("list", p)
  LSL_list <- vector("list", p)

  for(j in seq_len(p)){
    USL_list[[j]] <- funData(argvals = simData[[j]]@argvals, X = matrix(USL[[j]], nrow = 1))
    LSL_list[[j]] <- funData(argvals = simData[[j]]@argvals, X = matrix(LSL[[j]], nrow = 1)
    )
  }

  USL_fun <- multiFunData(USL_list)
  LSL_fun <- multiFunData(LSL_list)

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

  mfpca <- analisis$mfpca
  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)
  scores_mult <- analisis$scores_mult
  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 = M, ncol = 3)

  for(i in 1:M){
    cuantiles_proc[i, ] <- cuantiles(scores_mult[, i])
  }

  indices_cpk <- numeric(M)
  for(i in 1:M){
    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/M)
  return(list(lower = lower, upper = upper, indices_cpk = indices_cpk, MCnpk = MCnpk
  ))
}

### Función para graficar los scores ---------------------

graf_scores <- function(analisis, capacidad, M = 3){
  scores_mult <- analisis$scores_mult
  lower <- capacidad$lower
  upper <- capacidad$upper
  par(mfrow = c(1, M))
  for(m in 1:M){
    plot(scores_mult[,m], pch = 16, cex = 0.8,col = "black",
         xlab = "Observación", ylab = bquote(rho[.(m)]),
         main = paste("Score", m), ylim = c(min(min(scores_mult[,m]),lower[m]), max(max(scores_mult[,m]),upper[m])))
    abline(h = lower[m], col = "blue", lwd = 2, lty = 2)
    abline(h = upper[m], col = "red", lwd = 2, lty = 2)
  }
}

### Funcion para graficar los datos simulados con las especificaciones

graf_sim2 <- function(sim_obj, USL = NULL, LSL = NULL, ylim = NULL){
  simData <- sim_obj$simData
  p <- length(simData)
  X_list <- lapply(simData, function(fd) fd@X)
  if(is.null(ylim)){
    ymin_curvas <- min(sapply(X_list, min))
    ymax_curvas <- max(sapply(X_list, max))
    ymin <- ymin_curvas
    ymax <- ymax_curvas
    if(!is.null(LSL)){
      ymin <- min(ymin, min(unlist(LSL)))
    }
    if(!is.null(USL)){
      ymax <- max(ymax, max(unlist(USL)))
    }
    ylim <- c(ymin, ymax)
  }

  par(mfrow = c(1, p))

  for(j in seq_len(p)){
    t_grid <- simData[[j]]@argvals[[1]]
    Xj <- simData[[j]]@X

    matplot(t_grid, t(Xj), type = "l", lty = 1, col = "grey", xlab = "t", ylab = "", 
            main = paste("Componente", j), ylim = ylim)
    lines(t_grid, colMeans(Xj), lwd = 3, col = "black")
    if(!is.null(USL)){
      lines(t_grid, USL[[j]], col = "red", lwd = 3)
    }
    if(!is.null(LSL)){
      lines(t_grid, LSL[[j]], col = "blue", lwd = 3)
    }
  }
}

USL_sim2 <- list(USL1, USL2)
LSL_sim2 <- list(LSL1, LSL2)
Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim21 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "linear", N = 250, 
                    correlations = 0.1)


graf_sim2(sim_obj = sim21, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis21 <- analisis_mfpca2(sim_obj = sim21, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis21$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis21$UCL_T2))
abline(h = analisis21$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis21$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis21$UCL_Q))
abline(h = analisis21$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis21$reconstruction)

Code 
capacidad21 <- indices_capacidad_mfpca(sim_obj = sim21, analisis = analisis21,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad21$indices_cpk
[1] 1.8052313 0.1204218 0.7505375
Code 
capacidad21$MCnpk
[1] 0.5464328
Code 
graf_scores(analisis21, capacidad21, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim22 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "linear", N = 250, 
                    correlations = 0.4)


graf_sim2(sim_obj = sim22, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis22 <- analisis_mfpca2(sim_obj = sim22, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis22$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis22$UCL_T2))
abline(h = analisis22$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis22$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis22$UCL_Q))
abline(h = analisis22$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis22$reconstruction)

Code 
capacidad22 <- indices_capacidad_mfpca(sim_obj = sim22, analisis = analisis22,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad22$indices_cpk
[1] 1.7010183 0.3218521 0.2217070
Code 
capacidad22$MCnpk
[1] 0.495125
Code 
graf_scores(analisis22, capacidad22, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim23 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "linear", N = 250, 
                    correlations = 0.7)


graf_sim2(sim_obj = sim23, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis23 <- analisis_mfpca2(sim_obj = sim23, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis23$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis23$UCL_T2))
abline(h = analisis23$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis23$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis23$UCL_Q))
abline(h = analisis23$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis23$reconstruction)

Code 
capacidad23 <- indices_capacidad_mfpca(sim_obj = sim23, analisis = analisis23,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad23$indices_cpk
[1]  1.50447995  0.02095069 -0.01153180
Code 
capacidad23$MCnpk
[1] NaN
Code 
graf_scores(analisis23, capacidad23, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim24 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "linear", N = 250, 
                    correlations = 0.9)


graf_sim2(sim_obj = sim24, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis24 <- analisis_mfpca2(sim_obj = sim24, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis24$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis24$UCL_T2))
abline(h = analisis24$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis24$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis24$UCL_Q))
abline(h = analisis24$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis24$reconstruction)

Code 
capacidad24 <- indices_capacidad_mfpca(sim_obj = sim24, analisis = analisis24,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad24$indices_cpk
[1] 1.40880771 0.10856951 0.03470325
Code 
capacidad24$MCnpk
[1] 0.1744389
Code 
graf_scores(analisis24, capacidad24, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim25 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "exponential", N = 250, 
                    correlations = 0.1)


graf_sim2(sim_obj = sim25, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis25 <- analisis_mfpca2(sim_obj = sim25, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis25$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis25$UCL_T2))
abline(h = analisis25$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis25$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis25$UCL_Q))
abline(h = analisis25$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis25$reconstruction)

Code 
capacidad25 <- indices_capacidad_mfpca(sim_obj = sim25, analisis = analisis25,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad25$indices_cpk
[1]  1.72783699 -0.01309271  0.10883342
Code 
capacidad25$MCnpk
[1] NaN
Code 
graf_scores(analisis25, capacidad25, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim26 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "exponential", N = 250, 
                    correlations = 0.3)


graf_sim2(sim_obj = sim26, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis26 <- analisis_mfpca2(sim_obj = sim26, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis26$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis26$UCL_T2))
abline(h = analisis26$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis26$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis26$UCL_Q))
abline(h = analisis26$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis26$reconstruction)

Code 
capacidad26 <- indices_capacidad_mfpca(sim_obj = sim26, analisis = analisis26,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad26$indices_cpk
[1] 1.79309397 0.02712246 0.34650039
Code 
capacidad26$MCnpk
[1] 0.2563767
Code 
graf_scores(analisis26, capacidad26, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim27 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "exponential", N = 250, 
                    correlations = 0.8)


graf_sim2(sim_obj = sim27, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis27 <- analisis_mfpca2(sim_obj = sim27, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis27$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis27$UCL_T2))
abline(h = analisis27$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis27$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis27$UCL_Q))
abline(h = analisis27$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis27$reconstruction)

Code 
capacidad27 <- indices_capacidad_mfpca(sim_obj = sim27, analisis = analisis27,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad27$indices_cpk
[1] 1.42859451 0.16627285 0.07208267
Code 
capacidad27$MCnpk
[1] 0.2577431
Code 
graf_scores(analisis27, capacidad27, M=3)

Code 
t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim28 <- simMultiCorr(type = "split", argvals = argvals, M = 3, 
                    Ms= c(3,3), eFunType = "Fourier", 
                    eValType = "exponential", N = 250, 
                    correlations = 0.90)


graf_sim2(sim_obj = sim28, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis28 <- analisis_mfpca2(sim_obj = sim28, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis28$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis28$UCL_T2))
abline(h = analisis28$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis28$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis28$UCL_Q))
abline(h = analisis28$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis28$reconstruction)

Code 
capacidad28 <- indices_capacidad_mfpca(sim_obj = sim28, analisis = analisis28,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad28$indices_cpk
[1] 1.39074224 0.12825031 0.08157224
Code 
capacidad28$MCnpk
[1] 0.244127
Code 
graf_scores(analisis28, capacidad28, M=3)

Resultados

Code 
resultados_dependencia2 <- data.frame(
  
  Escenario = c(
    "Escenario 1",
    "Escenario 2",
    "Escenario 3",
    "Escenario 4",
    "Escenario 5",
    "Escenario 6",
    "Escenario 7",
    "Escenario 8"
  ),
  
  eValType = c(
    "linear",
    "linear",
    "linear",
    "linear",
    "exponential",
    "exponential",
    "exponential",
    "exponential"
  ),
  
  Correlación = c(0.1, 0.4, 0.7, 0.9, 0.1, 0.4, 0.7, 0.9),
  
  Cpk_1 = c(
    capacidad21$indices_cpk[1],
    capacidad22$indices_cpk[1],
    capacidad23$indices_cpk[1],
    capacidad24$indices_cpk[1],
    capacidad25$indices_cpk[1],
    capacidad26$indices_cpk[1],
    capacidad27$indices_cpk[1],
    capacidad28$indices_cpk[1]
    
  ),
  
  Cpk_2 = c(
    capacidad21$indices_cpk[2],
    capacidad22$indices_cpk[2],
    capacidad23$indices_cpk[2],
    capacidad24$indices_cpk[2],
    capacidad25$indices_cpk[2],
    capacidad26$indices_cpk[2],
    capacidad27$indices_cpk[2],
    capacidad28$indices_cpk[2]
  ),
  
  Cpk_3 = c(
    capacidad21$indices_cpk[3],
    capacidad22$indices_cpk[3],
    capacidad23$indices_cpk[3],
    capacidad24$indices_cpk[3],
    capacidad25$indices_cpk[3],
    capacidad26$indices_cpk[3],
    capacidad27$indices_cpk[3],
    capacidad28$indices_cpk[3]
  ),
  
  MCpk = c(
    capacidad21$MCnpk,
    capacidad22$MCnpk,
    capacidad23$MCnpk,
    capacidad24$MCnpk,
    capacidad25$MCnpk,
    capacidad26$MCnpk,
    capacidad27$MCnpk,
    capacidad28$MCnpk
  )
)

resultados_dependencia[,3:6] <- round(
  resultados_dependencia[,4:7],
  3
)

knitr::kable(
  resultados_dependencia2,
  caption = "Índices de capacidad para distintos escenarios de dependencia"
)
Índices de capacidad para distintos escenarios de dependencia
Escenario eValType Correlación Cpk_1 Cpk_2 Cpk_3 MCpk
Escenario 1 linear 0.1 1.805231 0.1204218 0.7505375 0.5464328
Escenario 2 linear 0.4 1.701018 0.3218521 0.2217070 0.4951250
Escenario 3 linear 0.7 1.504480 0.0209507 -0.0115318 NaN
Escenario 4 linear 0.9 1.408808 0.1085695 0.0347033 0.1744389
Escenario 5 exponential 0.1 1.727837 -0.0130927 0.1088334 NaN
Escenario 6 exponential 0.4 1.793094 0.0271225 0.3465004 0.2563767
Escenario 7 exponential 0.7 1.428594 0.1662728 0.0720827 0.2577431
Escenario 8 exponential 0.9 1.390742 0.1282503 0.0815722 0.2441270

Simulación 3

Funciones utilizadas

Code 
library(fda)
library(MASS)
library(funData)

sim_multi <- function(N = 250, argvals, M = 3, Sigma, basis_type = "fourier"){

  library(fda)
  library(MASS)
  library(funData)
  p <- length(argvals)
  if(nrow(Sigma) != M || ncol(Sigma) != M){
    stop("Sigma debe ser M x M")
  }

  global_grid <- seq(0, p, length.out = 1000)

  if(basis_type == "fourier"){
    basis_global <- create.fourier.basis(rangeval = c(0, p), nbasis = M)

  } else if(basis_type == "bspline"){

    basis_global <- create.bspline.basis(rangeval = c(0, p), nbasis = M)

  } else {
    stop("basis_type debe ser 'fourier' o 'bspline'")
  }

  Phi_global <- eval.basis(global_grid, basis_global)[,1:M]
  M_real <- ncol(Phi_global)
  if(M_real != M){
    warning(
      paste(
        "La base generó", M_real,
        "funciones. Ajustando automáticamente."
      )
    )
    if(nrow(Sigma) < M_real){
      stop("Sigma tiene menor dimensión que la base generada")
    }
    Sigma <- Sigma[1:M_real, 1:M_real]
    M <- M_real
  }
  
  Psi <- vector("list", p)
  for(j in 1:p){
    idx_j <- which(global_grid >= (j-1) &
                     global_grid <= j)
    Psi_j <- Phi_global[idx_j, , drop = FALSE]
    signs <- sample(c(-1,1), M, replace = TRUE)
    Psi_j <- sweep(Psi_j, 2, signs, "*")
    t_j <- argvals[[j]]
    Psi_interp <- matrix(0, nrow = length(t_j), ncol = M)
    local_grid <- seq(0,1,length.out=nrow(Psi_j))
    for(m in 1:M){
      Psi_interp[,m] <- approx(x = local_grid, y = Psi_j[,m], xout = t_j,
                               rule = 2)$y
    }

    for(m in 1:M){
      norma <- sqrt(sum(Psi_interp[,m]^2) * mean(diff(t_j)))
      Psi_interp[,m] <- Psi_interp[,m] / norma
    }
    Psi[[j]] <- Psi_interp
  }

  scores <- MASS::mvrnorm(n = N, mu = rep(0, M), Sigma = Sigma)

  X <- vector("list", p)
  for(j in 1:p){
    X[[j]] <- scores %*% t(Psi[[j]])
  }

  simData <- vector("list", p)
  for(j in 1:p){
    simData[[j]] <- funData(argvals = list(argvals[[j]]), X = X[[j]])
  }
  
  trueFuns <- vector("list", p)

  for(j in 1:p){
    trueFuns[[j]] <- funData(argvals = list(argvals[[j]]), X = t(Psi[[j]]))
  }

  trueVals <- eigen(Sigma)$values

  return(list(simData = multiFunData(simData), scores = scores, trueFuns = trueFuns,
              trueVals = trueVals, argvals = argvals, X = X))
}
2+1
[1] 3

Escenario 1

Code 
set.seed(NULL)

Sigma1 <- matrix(c(
  
  0.5, 0.0, 0.0,  0.08, 0.03, 0.01,
  0.0, 0.3, 0.0,  0.03, 0.07, 0.02,
  0.0, 0.0, 0.1,  0.01, 0.02, 0.04,
  
  0.08, 0.03, 0.01,  0.5, 0.0, 0.0,
  0.03, 0.07, 0.02,  0.0, 0.3, 0.0,
  0.01, 0.02, 0.04,  0.0, 0.0, 0.1
  
), nrow = 6, byrow = TRUE)

t1 <- seq(0, 1, length.out = 200)
t2 <- seq(0, 1, length.out = 200)
argvals <- list(t1, t2)

sim31 <- sim_multi(N=250, argvals = argvals, M=6, Sigma=Sigma1, basis_type = "fourier")
graf_sim2(sim_obj = sim31, USL = USL_sim2, LSL = LSL_sim2)

Code 
analisis31 <- analisis_mfpca2(sim_obj = sim31, alpha = 0.01, ncomp = 3, M = 3, R = 1000)

# Carta T2
plot(analisis31$T2, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = expression(T^2),
     main = expression("Carta " * T^2), ylim = c(0, analisis31$UCL_T2))
abline(h = analisis31$UCL_T2, col = "red", lwd = 2, lty = 2)

# Carta Q
plot(analisis31$Q, type = "p", pch = 16, cex = 0.7, xlab = "Observación", ylab = "Q", main = "Carta Q",
     ylim=c(0, analisis31$UCL_Q))
abline(h = analisis31$UCL_Q, col = "red", lwd = 2, lty = 2)

Code 
plot(analisis31$reconstruction)

Code 
capacidad31 <- indices_capacidad_mfpca2(sim_obj = sim31, analisis = analisis31,  
                                      USL = USL_sim2, LSL = LSL_sim2, M = 3)

capacidad31$indices_cpk
[1] 0.9384136 1.6324881 0.3899048
Code 
capacidad31$MCnpk
[1] 0.8421723
Code 
graf_scores(analisis31, capacidad31, M=3)

Code 
scores_corr(analisis31$scores_mult)

# scores univariados
scores_uni31 <- cbind(analisis31$pca_uni[[1]]$scores, analisis31$pca_uni[[2]]$scores)
corrplot(cor(scores_uni31), method = "color", addCoef.col = "black", tl.col = "black",
         tl.srt = 45)

Referencias

Happ, Clara, and Sonja Greven. 2018. “Multivariate Functional Principal Component Analysis for Data Observed on Different (Dimensional) Domains.” Journal of the American Statistical Association 113 (522): 649–59. https://doi.org/10.1080/01621459.2016.1273115.