Fluorescencia
La fluorescencia es la propiedad de algunos átomos y moléculas de absorber luz a una longitud de onda determinada (excitación) seguido por la emisión de luz de corta duración a una longitud de onda más larga.
La fluorescencia implica una fuente de luz externa para excitar la muestra a una longitud de onda en particular. Cuando se excita a la longitud de onda adecuada, la molécula pasa de un estado básico a otro excitado. A medida que la molécula vuelve al estado fundamental, se libera energía en forma de calor (pérdida de energía) y luz a una longitud de onda diferente de menor energía.
Recolección de los datos
Se recolectaron muestras de azúcar durante tres meses de operación en una planta azucarera ubicada en Escandinavia. Las muestras se tomaron cada ocho horas, lo que corresponde a tres mediciones diarias, para un total de 268 muestras a lo largo del periodo de estudio.
El azúcar se disolvió en agua en una proporción de 2.25 g por 15 ml, y la solución resultante se analizó utilizando un espectrofluorómetro PE LS50B.
Cada muestra fue excitada a siete longitudes de onda (230, 240, 255, 290, 305, 325 y 340 nm). Para cada una de estas excitaciones, el equipo registró la intensidad de emisión en 571 longitudes de onda. Los espectros de emisión se midieron en el rango de 275 a 560 nm, con intervalos de 0,5 nm.
Para el suavizado de las curvas se empleará suaviado splines, utilizando bases B-splines. Esta elección se justifica a partir de la visualización de los datos, donde se puede observar que los espectros de emisión no presentan un comportamiento periódico. En este contexto, los B-splines son una base adecuada.
#| code-fold: false
#| echo: true
file <- read.mat("data.mat")
data <- file$X
t <- file$EmAx
wvl <- array(NA, dim = c(571, 268, 7))
for (k in 1:7) {
start <- 571 * (k - 1) + 1
end <- 571 * k
wvl[, , 8 - k] <- t(data[,start:end])
}
nms <- c("340nm", "325nm", "305nm", "290nm", "255nm", "240nm", "230nm")
nms <- rev(nms)
dimnames(wvl) <- list(NULL, NULL, nms)Basis <- create.bspline.basis(rangeval = c(275, 560),
norder = 6, nbasis = 40)
fit <- list()
grid <- seq(min(t), max(t), length = 1000)
fits <- array(dim = c(length(grid), dim(wvl)[2:3]))
for (k in 1:7) {
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = Lambdas[k])
fit[[k]] <- smooth.basis(t, wvl[, , k], fdParobj)
fdobj <- fit[[k]]$fd
fits[ , , k] <- eval.fd(evalarg = grid, fdobj = fdobj)
}
Para poder realizar FPCA, el operador de covarianza debe ser un operador de Hilbert–Schmidt, simétrico y definido positivo. Para que lo primero se cumpla, el proceso aleatorio del cual se obtienen las realizaciones debe ser cuadrado–integrable, \(\mathrm{E}[\lVert \chi \rVert^2]<\infty\). Las otras dos propiedades se tienen por la simetría y positividad del operador de covarianza \(c(s,t) = \operatorname{Cov}(\chi(s), \chi(t))\).
Otro supuesto que se hace es que la media es constante. Sean \(\chi_1\),\(\chi_2\)\(\cdots\)\(\chi_N\) observaciones funcionales univariadas en \(L^2\).
Para esto, se juzga el sistema de hipótesis \[ H_0: \mathrm{E}[\chi_1] = \cdots = \mathrm{E}[\chi_N] \] Cada observación funcional puede representarse mediante el modelo:
\[ \begin{equation} \chi_i(t) = \mu(t) + Y_i(t), \quad E[Y_i(t)] = 0 \end{equation} \] Donde \(\mu(t)\) es la media poblacional.
Estadístico de Prueba
Para explicar la idea del procedimiento de prueba, se define, de acuerdo con (Horváth y Kokoszka 2012) y (Berkes et al. 2009)
\[ \hat{\mu}_k(t) = \frac{1}{k} \sum_{1 \leq i \leq k} \chi_i(t), \qquad \check{\mu}_k(t) = \frac{1}{N-k} \sum_{k < i \leq N} \chi_i(t). \]
Si la media es constante, la diferencia
\[ \Delta_k(t) = \hat{\mu}_k(t) - \check{\mu}_k(t) \]
es pequeña para todo \(1 \leq k < N\) y para todo \(t \in [0,1]\).
Sin embargo, \(\Delta_k(t)\) puede volverse grande debido a la variabilidad si \(k\) está cerca de \(1\) o de \(N\). Por esta razón, es usual trabajar con
Si \(k\) es cercano a 1, \(\hat{\mu}_k\) se calcula con pocos valores, mientras que \(\tilde{\mu}_c¿k\) incorpora la mayor parte de la información. En cambio, si \(k\) es cercano a \(N\), ocurre lo contrario. En cualquiera de estos dos casos, \(\hat{\mu}_k\) y \(\tilde{\mu}_k\) pueden diferir considerablemente, lo que podría llevar a que la prueba detecte un punto de cambio incluso cuando en realidad no existe.
Para solucionar este problema, se puede ponderar la diferencia entre estos valores mediante un peso parabólico en función de \(k\). \[ P_k(t) = \sum_{1 \leq i \leq k} \chi_i(t) - \frac{k}{N} \sum_{1 \leq i \leq N} \chi_i(t) = \frac{k(N-k)}{N} \left[ \hat{\mu}_k(t) - \check{\mu}_k(t) \right]. \]
Ya que el valor de \(P_k(t)\) no cambia si las funciones \(\chi_i(t)\) son reemplazadas por \(\chi_i(t)-\chi_N(t)\). Al definir \(k=[Nx]\), con \(x \in (0,1)\), se obtiene
\[ \begin{equation} \int \left\{ \sum_{1 \leq i \leq Nx} \chi_i(t) - \frac{[Nx]}{N} \sum_{1 \leq i \leq N} \chi_i(t) \right\} \hat{v}_{\ell}(t)\,dt = \sum_{1 \leq i \leq Nx} \hat{\xi}_{\ell,i} - \frac{[Nx]}{N} \sum_{1 \leq i \leq N} \hat{\xi}_{\ell,i}. \end{equation} \]
Lo que muestra que los scores pueden utilizarse para probar la constancia de la función de media.
El siguiente teorema puede emplearse para derivar varios estadísticos de prueba.
\[ \begin{equation} T_N(x) = \frac{1}{N} \sum_{\ell=1}^{d} \hat{\lambda}_{\ell}^{-1} \left( \sum_{1 \leq i \leq Nx} \hat{\xi}_{\ell,i} - x \sum_{1 \leq i \leq N} \hat{\xi}_{\ell,i} \right)^2 \end{equation} \]
Sean \(B_1(\cdot), \ldots, B_d(\cdot)\) puentes brownianos estándar independientes. Bajo \(H_0\):
\[ T_N(x) \overset{d}{\longrightarrow} \sum_{1 \leq \ell \leq d} B_{\ell}^{2}(x), \qquad (0 \leq x \leq 1), \]
El estadístico de prueba tipo Cramér-von Mises, denotado como \(S_{N,d}\), se define para \(d\) componentes principales como:
\[ \begin{equation} S_{N,d} = \frac{1}{N^2} \sum_{\ell=1}^{d} \frac{1}{\hat{\lambda}_\ell} \sum_{k=1}^{N} \left( \sum_{i=1}^k \hat{\xi}_{\ell,i} - \frac{k}{N} \sum_{i=1}^N \hat{\xi}_{\ell,i} \right)^2 \end{equation} \]
Bajo la hipótesis nula, el estadístico \(S_{N,d}\) converge en distribución a la suma de cuadrados de puentes brownianos independientes:
\[ \begin{equation} S_{N,d} \xrightarrow{d} \int_0^1 \sum_{\ell=1}^d B_\ell^2(x) dx \end{equation} \]
La decisión estadística se toma comparando el valor calculado de \(S_{N,d}\) con los valores críticos. Se rechaza \(H_0\) con un nivel de significancia \(\alpha\) si: \[ \begin{equation} S_{N,d} > C_d(\alpha) \end{equation} \] Teniendo en cuenta que \(x \in [0, 1]\), el valor \(x^*\) que maximiza \(T_N(x)\), multiplicado por \(N\) corresponde a un estimador consistente del punto de cambio.
Esto se realiza mediante la función fchange() del paquete fChange, donde es posible aplicar la metodología propuesta por (Berkes et al. 2009). Esta metodología
X.dfts <- dfts(fits[ , , 2])
N <- ncol(X.dfts$data)
ChangeDetectionSeg <- fchange(X.dfts,
method = "projmean",
statistic = "Tn",
type = "segmentation",
critical = "simulation")| Punto de Cambio | p-valor |
|---|---|
| 3 | 0.038 |
| 7 | 0.000 |
| 13 | 0.043 |
| 16 | 0.005 |
| 112 | 0.003 |
| 134 | 0.000 |
| 188 | 0.000 |
| 194 | 0.007 |
| 204 | 0.002 |
| 216 | 0.038 |
| 246 | 0.005 |
| 256 | 0.014 |
Los resultados proporcionados en Tabla 1 muestran evidencia para rechazar la hipótesis nula de media funcional constante, \(H_0: \mathrm{E}[\chi_1] = \cdots = \mathrm{E}[\chi_N]\), a favor de la existencia de cambios en la media de las curvas espectrométricas.
| Segmento | n | Rango |
|---|---|---|
| 1 | 3 | 1 a 3 |
| 2 | 4 | 4 a 7 |
| 3 | 6 | 8 a 13 |
| 4 | 3 | 14 a 16 |
| 5 | 96 | 17 a 112 |
| 6 | 22 | 113 a 134 |
| 7 | 54 | 135 a 188 |
| Segmento | n | Rango | |
|---|---|---|---|
| 8 | 8 | 6 | 189 a 194 |
| 9 | 9 | 10 | 195 a 204 |
| 10 | 10 | 12 | 205 a 216 |
| 11 | 11 | 30 | 217 a 246 |
| 12 | 12 | 10 | 247 a 256 |
| 13 | 13 | 12 | 257 a 268 |
Por lo tanto, a cada observación se le resta la media estimada correspondiente al segmento al que pertenece, con el fin de obtener las curvas centradas.
A continuación se realiza el análisis de componentes principales para el segmento 5 y el segmento 7.
X.G1 <- cp[[idxSeg[1]]]$fd_original
pcaObj.G1 <- pca.fd(X.G1, nharm = 3, centerfns = TRUE)
eigenfn.G1 <- eval.fd(grid, pcaObj.G1$harmonics)
eigenval.G1 <- pcaObj.G1$values[1:3]X.G2 <- cp[[idxSeg[2]]]$fd_original
pcaObj.G2 <- pca.fd(X.G2, nharm = 3, centerfns = TRUE)
eigenfn.G2 <- eval.fd(grid, pcaObj.G2$harmonics)
eigenval.G2 <- pcaObj.G2$values[1:3]Luego de haber centrado los procesos, se encuentran las estimaciones de las funciones propias asociadas al operador de covarianza del proceso funcional para cada grupo.
[1] "done"
Con la función plot.pca.fd es posible interpretarlas de una forma más sencilla, pues se muestra cómo cada componente principal afecta a la media del proceso. Esto se muestra en la Figura 4, Figura 5 y Figura 6.
[1] "done"
En el caso de la primera función propia, esta está mostrando un aumento constante en la media. Debido a su alto porcentaje de variabilidad explicada, se podría decir que la mayor fuente de variación es la cantidad de luz que se emite (magnitud). Además, se puede ver que hay más varibilidad cuando las muestras se someten a longitudes de onda entre los 350nm y 500nm aproximadamente.
La segunda función propia muestra un contraste entre la exposición a las longitudes de onda menores a 380nm aproximadamente y la exposición a las longitudes de onda mayores a 380nm, mostrando una variación mayor en los intervalos [310, 370] y [400, 470] nm.
Por último, la tercera función propia tienen un patrón cuadrático y presenta mayor variabilidad hacia las longitudes de onda medias.
El porcentaje de varianza se muestra a través de un scree plot en la Figura 11. De esta forma, es posible ver que para los dos segmentos, las dos primeras componentes principales ya explican más del 90% de variabilidad total.
La descomposición de Karhunen-Loève consiste en expresar cada observación funcional como \[ X(t) = \mu(t) + \sum_{k=1}^\infty \xi_k \nu_k (t) \] donde \(\nu_k\) corresponde a la \(k\)–ésima función propia y \(\xi_k\) es el \(k\)–ésimo score, \(\xi_k=\langle X, \nu_k \rangle\).
En este caso, la base óptima corresponde al conjunto de las dos primeras funciones propias pues ellas explican más del 90% de la variabilidad, aproximadamente. Por lo tanto, para cada observación funcional, se tiene su aproximación \[ X_{i,g}(t) \approx \hat \mu_g(t) + \sum_{k=1}^{K_g}\hat {\xi}_{ikg} \hat \nu _ {kg}(t) \] con \(i = 1,\ldots,268\), \(g=1,2\) y \(K_g=2\).
A continuación, se expresan las tres primeras observaciones muestrales con la descomposición de Karhunen-Loève.
K1 <- 2
Scores.G1 <- matrix(pcaObj.G1$scores[,1:K1], nrow = nrow(X.G1), ncol = K1)
eigenfns.G1 <- eval.fd(evalarg = grid, fdobj = pcaObj.G1$harmonics)[,1:K1]
mu.G1 <- eval.fd(evalarg = grid, fdobj = pcaObj.G1$meanfd)
X.approx.G1 <- Scores.G1 %*% t(eigenfns.G1)
X.approx.G1 <- sweep(X.approx.G1, 2, mu.G1, "+") 
K2 <- 2
Scores.G2 <- matrix(pcaObj.G2$scores[,1:K2], nrow = nrow(X.G2), ncol = K2)
eigenfns.G2 <- eval.fd(evalarg = grid, fdobj = pcaObj.G2$harmonics)[,1:K2]
mu.G2 <- eval.fd(evalarg = grid, fdobj = pcaObj.G2$meanfd)
X.approx.G2 <- Scores.G2 %*% t(eigenfns.G2)
X.approx.G2 <- sweep(X.approx.G2, 2, mu.G2, "+") 
A continuación, en la Tabla 2 se muestra el error cuadrático medio de predicción para cada observación.
| MPSE | |
|---|---|
| Segmento 5 | 389.2850 |
| Segmento 7 | 788.0029 |
En el caso del primer grupo es posible ver que los scores tienen una distribución asimétrica y también hay presencia de algunos outliers.
A continuación, se muestra el gráfico de dispersión de las dos primeras componentes principales, que explican más del 90% de la variabilidad, aproximadamente.
Es posible ver que la mayoría de puntos se concentran alrededor de los valores negativos cercanos al cero para el primer score y alrededor del cero para el segundo score. Asimismo, debido al gran porcentaje de variabilidad que explica la primera componente se puede ver que la dispersión es mayor a lo largo de su eje.
A partir de la gráfica, se observa que la función de covarianza estimada alcanza sus valores más altos en las longitudes de onda intermedias. Esto indica que la variabilidad de las curvas es mayor en la región central.
Este comportamiento es coherente con lo observado en la primera función propia, la cual presenta sus valores más altos precisamente en dichas longitudes de onda, lo que sugiere que la mayor parte de la variabilidad está concentrada en esa zona y corresponde a un patrón global de cambio en el nivel de las curvas.
La segunda y tercera funciones propias muestran cambios de signo y formas más oscilatorias, lo que refleja variaciones más complejas en la forma de las curvas. Estas funciones permiten capturar patrones de contraste y variaciones locales que no son explicados por la primera componente principal, aportando así mayor flexibilidad en la representación de la variabilidad funcional.
Por el teorema de Mercer se tiene que el kernel \(\psi(s,t)\) usado en el operador integral \(\Psi (x)(t)\) se puede representar en términos de los valores propios y funciones propias del operador. En este caso, se tiene el operador de covarianza \(\mathcal{C}_\chi(x)(t))\) con la función de covarianza como Kernel. Luego,
\[c(s,t) = \sum_{j=1}^\infty \lambda_j \phi(s) \phi(t)\]
A continuación se puede ver la función de covarianza junto con su aproximación y sus residuales.
Una carta de control es una herramienta que permite monitorear si un proceso se mantiene bajo control o si presenta cambios anómalos a lo largo del tiempo.
Para construir una carta de control a partir del análisis de componentes principales:
Para el caso funcional, sea \(\chi\) una función aleatoria definida en \(L^2\), la cual puede representarse mediante la descomposición de Karhunen–Loève como:
\[\begin{equation} \chi(t) = \mu(t) + \sum_{k=1}^{\infty} \xi_k \, \phi_k(t) \end{equation}\]
Donde:
Para realizar FPCA, se utiliza una aproximación truncada de esta expansión. Dado que el operador de covarianza posee infinitos valores y funciones propias, se selecciona un número finito \(K\) de componentes principales, de manera que se capture una proporción suficientemente alta de la varianza total.
Así, se obtiene la aproximación:
\[\begin{equation} \chi(t) \approx \mu(t) + \sum_{k=1}^{K} \xi_k \, \phi_k(t) \end{equation}\]
La carta \(T^2\) se contruye sobre los scores de la siguiente manera:
\[\boldsymbol{\xi} = \left( \langle \chi - \mu, \phi_1 \rangle,\; \langle \chi - \mu, \phi_2 \rangle,\; \ldots,\; \langle \chi - \mu, \phi_K \rangle \right)^{'}\] Donde
\[\begin{align} \operatorname{Cov}(\xi_i, \xi_j) &= \operatorname{Cov}\big(\langle \chi - \mu, \phi_i \rangle,\; \langle \chi - \mu, \phi_j \rangle \big) \\ &= 0 \quad \text{si } i \neq j \\[6pt] \operatorname{Var}(\xi_i) &= \operatorname{Cov}(\xi_i, \xi_i) \\ &= \operatorname{Cov}\big(\langle \chi - \mu, \phi_i \rangle,\; \langle \chi - \mu, \phi_i \rangle \big) \\ &= \lambda_i \end{align}\]
Por lo tanto \(\Sigma_{\boldsymbol{\xi}} = \operatorname{diag}(\lambda_1, \lambda_2, \ldots, \lambda_K)\)
Sea \(x_1,x_2,..,x_n\) realizaciones de \(\chi\) tal que:
\[\boldsymbol{\xi}_j = \left( \langle x_j - \bar{x}, \phi_1 \rangle,\; \langle x_j - \bar{x}, \phi_2 \rangle,\; \ldots,\; \langle x_j - \bar{x}, \phi_K \rangle \right)^{'}.\]
Entonces el estadístico $ T_j^2$ se define como:
\[\begin{equation} T_j^2 = \boldsymbol{\xi}_j^{\top}\,\Sigma_{\boldsymbol{\xi}}^{-1}\,\boldsymbol{\xi}_j = \sum_{k=1}^{K} \frac{\xi_{jk}^2}{\lambda_k}, \end{equation}\]
El estadístico \(Q_j\) se define como:
\[ \begin{equation} Q_j = \| \varepsilon_j \|^2 =\| x_j -\widehat{x_j}\|^2= \left\| x_j - \bar{x} - \sum_{k=1}^{K} \xi_{jk}\,\phi_k \right\|^2, \end{equation} \]
En la Fase I se realiza el FPCA a partir de una muestra en control.
El nivel de significancia conjunto para las cartas de control \(T^2\) y \(Q\) se fija como \(\alpha\). Dado que se utilizan dos estadísticas de monitoreo, el nivel de significancia individual para cada carta se define como:
\[ \alpha' = 1 - (1 - \alpha)^{1/p}, \]
donde \(p = 2\) corresponde al número de cartas de control. De esta manera, se garantiza que el nivel de significancia global sea \(\alpha\). Nosotros fijaremos \(\alpha = 0.05\) como nivel de significancia conjunto.
Por otra parte, dado que la distribución exacta de los estadísticos \(T^2\) y \(Q\) no es conocida, los límites de control se estiman mediante un procedimiento de bootstrap. En particular, se considera el máximo de estos estadísticos en cada réplica bootstrap, y a partir de esta distribución empírica se calculan los cuantiles correspondientes de \(T^2\) y \(Q\), los cuales se utilizan para definir los limites.
\[ \text{LC}_{T^2} = q_{1-\alpha'}\left( {T_1^2}, \ldots, T_n^2 \right), \]
\[ \text{LC}_{Q} = q_{1-\alpha'}\left( Q_1, \ldots, Q_n \right), \]
Este enfoque permite evitar supuestos distribucionales y proporciona límites de control adaptados a la estructura de los datos.
Al final de la Fase I, se debe quedar con el conjunto de datos que represente un proceso en control.
Para una nueva observación \(x_{\text{new}}\), se calculan:
\[ T^2 = \sum_{k=1}^{K} \frac{\xi_k^2}{\lambda_k}, \qquad Q = \left\| x_{\text{new}} - \bar{x} - \sum_{k=1}^{K} \xi_k\,\phi_k \right\|^2. \]
Se considera fuera de control si:
\[ T^2 > \text{UCL}_{T^2} \quad \text{o} \quad Q > \text{UCL}_{Q}. \]
Para la construcción de las cartas de control, se divide el conjunto de datos en dos subconjuntos: el \(80\%\) de las observaciones se utiliza en la Fase I, mientras que el \(20\%\) restante se reserva para la Fase II.
Una vez identificados los puntos de cambio, la construcción de las cartas de control se realizará mediante la partición del conjunto de datos, agrupando aquellas observaciones que comparten una misma estructura de media.
T2_Q.CL <- function(data, argvals, Lambda, Basis, ncomp, alpha, B, seed){
N <- dim(data)[2]
all_T2 <- vector("list", B)
all_Q <- vector("list", B)
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = Lambda)
dt <- diff(argvals)
set.seed(seed)
for (b in seq_len(B)) {
# Remuestreo
idx <- sample(1:N, size = N, replace = TRUE)
samples <- data[, idx, drop = FALSE]
# Suavizado y PCA
smoothlist <- smooth.basis(argvals, samples, fdParobj)
fdObj_c <- center.fd(smoothlist$fd)
pcaObj <- pca.fd(fdObj_c, nharm = ncomp, centerfns = FALSE)
# Extraer scores y valores propios
Scores <- pcaObj$scores
Eigenval <- pcaObj$values[1:ncomp]
# T2
all_T2[[b]] <- rowSums((Scores %*% diag(1/Eigenval)) * Scores)
# Q
Xfd <- eval.fd(argvals, fdObj_c)
Phi <- eval.fd(argvals, pcaObj$harmonics)
Xhat <- Phi %*% t(Scores)
all_Q[[b]] <- apply((Xfd - Xhat)^2, 2, function(v){
sum((v[-1] + v[-length(v)]) / 2 * dt)
})
}
# Límites de control
CL_T2 <- quantile(unlist(all_T2), 1 - alpha)
CL_Q <- quantile(unlist(all_Q), 1 - alpha)
# Modelo con datos originales
smoothlist_f <- smooth.basis(argvals, data, fdParobj)
fdObj_c_f <- center.fd(smoothlist_f$fd)
pcaObj_f <- pca.fd(fdObj_c_f, nharm = ncomp, centerfns = FALSE)
Scores_f <- pcaObj_f$scores
Eigenval_f <- pcaObj_f$values[1:ncomp]
# T2 final
T2_final <- rowSums((Scores_f %*% diag(1/Eigenval_f)) * Scores_f)
# Q final
Xfd_f <- eval.fd(argvals, fdObj_c_f)
Phi_f <- eval.fd(argvals, pcaObj_f$harmonics)
Xhat_f <- Phi_f %*% t(Scores_f)
Q_final <- apply((Xfd_f - Xhat_f)^2, 2, function(v){
sum((v[-1] + v[-length(v)]) / 2 * dt)
})
return(list(
CL_T2 = as.numeric(CL_T2),
CL_Q = as.numeric(CL_Q),
T2 = T2_final,
Q = Q_final
))
}
clean_phaseI_T2Q <- function(data, p, argvals, Lambdas, Basis, ncomp = 2,
alpha = 0.025, B = 500, seed = 123,r,cp){
set.seed(seed)
data_variable <- data[ , , p][,cp]
n_total <- dim(data_variable)[2]
# Selección de Fase I
idx_phase1 <- 1:floor(r * n_total)
data_phase1 <- data_variable[, idx_phase1, drop = FALSE]
indices <- idx_phase1
removed_indices <- c()
Lambda <- Lambdas[p]
repeat {
res <- T2_Q.CL(
data = data_phase1,
argvals = argvals,
Lambda = Lambda,
Basis = Basis,
ncomp = ncomp,
alpha = alpha,
B = B,
seed = seed
)
# Identificar fuera de control
outliers <- which(res$T2 > res$CL_T2 | res$Q > res$CL_Q)
if(length(outliers) == 0) break
# Guardar cuáles estamos quitando
removed_indices <- c(removed_indices, indices[outliers])
# Actualizar datos y los índices que quedan
data_phase1 <- data_phase1[, -outliers, drop = FALSE]
indices <- indices[-outliers]
if(length(indices) <= ncomp) break
}
return(list(
cleaned_data = data_phase1,
final_result = res,
removed_indices = removed_indices,
n_removed = length(removed_indices),
indices = indices
))
}Se construirá la carta de control para el proceso en el cual la media se mantiene constante y que, además, contiene el mayor número de observaciones. En este caso, corresponde al segmento 5.
La división de los datos se realizará utilizando el 80% de las observaciones para la Fase I y el 20% restante para la Fase II. Finalmente, los límites de control serán estimados mediante un procedimiento de bootstrap.
resT2Q_2=clean_phaseI_T2Q (data=wvl, p=2, argvals=t, Lambdas=Lambdas, Basis=Basis,B = 500,r=0.8,cp=17:112)Del segmento 5, conformado por 96 observaciones que presentan la misma media, se seleccionó el 80%, correspondiente a 76 observaciones, para la Fase I. Dentro de este conjunto, 41 observaciones se utilizaron para la estimación de los límites de control, mientras que las 20 observaciones restantes se destinaron a la Fase II, con el propósito de realizar el monitoreo del proceso.
A partir de la Fase I, los límites de control obtenidos fueron 5.866352 para el estadístico \(T^2\) y 173.5451 para el estadístico \(Q\).
dt <- diff(t)
# Scores y valores propios segmento 5
Scores_g3 <- pcaObj.G1$scores[,1:2]
Eigenval_g3 <- pcaObj.G1$values[1:2]
# Q segmento 5
Xfd_g3 <- eval.fd(t, cp$Grupo_5$fd_centrado)
Phi_g3 <- eval.fd(t, pcaObj.G1$harmonics)[,1:2]
Xhat_g3 <- Phi_g3 %*% t(Scores_g3)
indices_excluir <- c(
resT2Q_2[["indices"]],
resT2Q_2[["removed_indices"]]
)
# Índices válidos (los que NO están en el vector anterior)
indices_validos <- setdiff(
1:nrow(Scores_g3),
indices_excluir
)
# T2 solo para los índices válidos
T2_final_3 <- rowSums(
(Scores_g3[indices_validos, ] %*% diag(1 / Eigenval_g3)) *
Scores_g3[indices_validos, ]
)
# Q segmento 5 solo para los índices válidos
Q_final_3 <- apply(
(Xfd_g3[, indices_validos] - Xhat_g3[, indices_validos])^2,
2,
function(v) {
sum((v[-1] + v[-length(v)]) / 2 * dt)
}
)
CL_T2 <- resT2Q_2$final_result$CL_T2
CL_Q <- resT2Q_2$final_result$CL_Q
df_g3 <- data.frame(
obs = indices_validos + 16,
T2 = T2_final_3,
Q = Q_final_3
)
df_long_g3 <- df_g3 %>%
pivot_longer(
cols = c(T2, Q),
names_to = "stat",
values_to = "value"
)
limits <- data.frame(
stat = c("T2", "Q"),
CL = c(CL_T2, CL_Q)
)
palette <- c(
"T2" = "#1F77B4",
"Q" = "#FF7F0E"
)
limit_color <- "#D62728"
ggplot(df_long_g3, aes(x = obs, y = value, color = stat)) +
geom_line(linewidth = 0.9) +
geom_point(size = 1.6, alpha = 0.85) +
geom_hline(
data = limits,
aes(yintercept = CL),
color = limit_color,
linetype = "dashed",
linewidth = 1
) +
facet_wrap(~stat, ncol = 1, scales = "free_y") +
scale_color_manual(values = palette) +
labs(
x = "Índice de observación original",
y = "Valor del estadístico",
color = "Estadístico"
) +
theme_minimal(base_size = 13) +
theme(
legend.position = "none",
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", hjust = 0.5),
panel.grid.minor = element_blank()
)
| Observación real | Valor T² | Valor Q | Fuera por T² | Fuera por Q |
|---|---|---|---|---|
| 93 | 2.883 | 790.247 | No | Sí |
| 94 | 2.240 | 174.686 | No | Sí |
| 95 | 0.211 | 1780.197 | No | Sí |
| 96 | 1.723 | 1189.014 | No | Sí |
| 99 | 0.400 | 453.050 | No | Sí |
| 100 | 0.098 | 947.013 | No | Sí |
| 101 | 0.264 | 733.997 | No | Sí |
| 102 | 0.992 | 617.703 | No | Sí |
| 103 | 0.496 | 434.745 | No | Sí |
| 105 | 1.380 | 521.318 | No | Sí |
| 106 | 0.979 | 419.474 | No | Sí |
| 107 | 0.111 | 1324.682 | No | Sí |
| 108 | 2.029 | 521.007 | No | Sí |
| 109 | 6.747 | 2151.488 | Sí | Sí |
| 110 | 1.287 | 276.988 | No | Sí |
Sea \(\chi\) una función aleatoria definida en \(L^2\), asumimos \[ \mu_\chi(t) = E[\chi(t)], \] y \(\mathcal{C}_X\) el operador de covarianza
\[ \mathcal{C}_\chi(\eta) = E\big[(X - \mu_X) \otimes (X - \mu_X)\big] \]
Sean \(\lambda_1, \lambda_2, \ldots\) los valores propios del operador y \(\psi_1, \psi_2, \ldots\) las funciones propias, que forman una base ortonormal de \(L^2\).
Usando la expansión de Karhunen–Loève de la función aleatoria \(X\): \[ \chi = \mu_\chi + \sum_{k=1}^{\infty} \xi_k \psi_k\] donde
\[ \xi_k = \langle \chi - \mu_X, \psi_k \rangle, \] los cuales son los scores de componentes principales, variables aleatorias no correlacionadas con media cero y varianza \(\lambda_k\).
El operador de covarianza se expresa como: \[ \mathcal{C}_\chi(x,y) = \sum_{j=1}^{\infty} \lambda_j \psi_j(x)\psi_j(y), \] usando el teorema de Mercer.
Para una función \(y\), el operador inverso se escribe como: \[ \mathcal{C}_\chi^{-1}(y) = \sum_{k=1}^{\infty} \frac{1}{\lambda_k} \langle y, \psi_k \rangle \psi_k. \]
Se define la distancia funcional de Mahalanobis entre \(\chi\) y \(\mu_\chi\) como: \[ d_{FM}^2 = \left\langle \mathcal{C}_\chi^{-1/2}(\chi - \mu_\chi), \mathcal{C}_\chi^{-1/2}(\chi - \mu_\chi) \right\rangle. \]
Expandiendo: \[ \mathcal{C}_\chi^{-1/2}(\chi - \mu_\chi) = \sum_{k=1}^{\infty} \frac{1}{\sqrt{\lambda_k}} \xi_k \psi_k, \] por lo que
\[ d_{FM}^2 = \sum_{k=1}^{\infty} \frac{\xi_k^2}{\lambda_k}. \]
Equivalentemente,
\[ d_{FM} = \left( \sum_{k=1}^{K} \frac{\xi_k^2}{\lambda_k} \right)^{1/2}. \]
Esta versión funcional de la distancia de Mahalanobis es una semi-distancia, ya que debido al truncamiento del operador de covarianza puede ocurrir que elementos distintos tengan distancia cero.
Para la construcción de las cartas de control se seguirá un procedimiento análogo al utilizado previamente en el contexto de FPCA, dividiendo el análisis en dos etapas: una fase de calibración (Fase I) y una fase de monitoreo (Fase II).
En la fase I se seleccionará un subconjunto de observaciones que se asume proviene de un proceso bajo control, es decir, datos que comparten la misma estructura de media. En este caso, se considerará nuevamente que el periodo estable corresponde a las observaciones finales, por lo que se tomará el 80% de los datos más recientes para conformar la Fase I, mientras que el 20% restante se destinará a la Fase II.
Esta partición permitirá estimar los límites de control a partir de un conjunto de referencia estable y, posteriormente, evaluar el comportamiento de las observaciones restantes. Manteniendo el mismo enfoque basado en bootstrap y en la distancia funcional de Mahalanobis.
DFM.CL <- function(data, argvals, Lambda, Basis, ncomp, alpha, B, seed){
N <- dim(data)[2]
all_DFM <- vector("list", B)
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = Lambda)
set.seed(seed)
for (b in seq_len(B)) {
# Remuestreo bootstrap
idx <- sample(1:N, size = N, replace = TRUE)
samples <- data[, idx, drop = FALSE]
# Suavizado y FPCA
smoothlist <- smooth.basis(argvals, samples, fdParobj)
fdObj_c <- center.fd(smoothlist$fd)
pcaObj <- pca.fd(fdObj_c, nharm = ncomp, centerfns = FALSE)
# Scores y valores propios
Scores <- pcaObj$scores
Eigenval <- pcaObj$values[1:ncomp]
# Distancia funcional de Mahalanobis
all_DFM[[b]] <- sqrt(
rowSums(
(Scores^2) / matrix(
Eigenval,
nrow = nrow(Scores),
ncol = length(Eigenval),
byrow = TRUE
)
)
)
}
CL_DFM <- quantile(unlist(all_DFM), 1 - alpha)
# Modelo final con los datos originales
smoothlist_f <- smooth.basis(argvals, data, fdParobj)
fdObj_c_f <- center.fd(smoothlist_f$fd)
pcaObj_f <- pca.fd(fdObj_c_f, nharm = ncomp, centerfns = FALSE)
Scores_f <- pcaObj_f$scores
Eigenval_f <- pcaObj_f$values[1:ncomp]
# DFM final
DFM_final <- sqrt(
rowSums(
(Scores_f^2) / matrix(
Eigenval_f,
nrow = nrow(Scores_f),
ncol = length(Eigenval_f),
byrow = TRUE
)
)
)
return(list(
CL_DFM = as.numeric(CL_DFM),
DFM = DFM_final
))
}
clean_phaseI_DFM <- function(data, p, argvals,Lambdas,
Basis, ncomp = 2, alpha = 0.025, B = 500, seed = 123,
r, cp
) {
set.seed(seed)
data_variable <- data[, , p][, cp]
n_total <- dim(data_variable)[2]
# Selección aleatoria de Fase I
idx_phase1 <- 1:floor(r * n_total)
data_phase1 <- data_variable[, idx_phase1, drop = FALSE]
# Guardar índices originales
indices <- idx_phase1
removed_indices <- c()
Lambda <- Lambdas[p]
repeat {
# Estimar límite de control y DFM
res <- DFM.CL(
data = data_phase1,
argvals = argvals,
Lambda = Lambda,
Basis = Basis,
ncomp = ncomp,
alpha = alpha,
B = B,
seed = seed
)
# Identificar observaciones fuera de control
outliers <- which(res$DFM > res$CL_DFM)
# Si no hay outliers, detener
if (length(outliers) == 0) break
# Guardar índices eliminados
removed_indices <- c(
removed_indices,
indices[outliers]
)
# Eliminar observaciones fuera de control
data_phase1 <- data_phase1[, -outliers, drop = FALSE]
indices <- indices[-outliers]
# Evitar que queden muy pocas observaciones
if (length(indices) <= ncomp) break
}
return(list(
cleaned_data = data_phase1,
final_result = res,
removed_indices = removed_indices,
n_removed = length(removed_indices),
indices = indices
))
}
CL_DFM=clean_phaseI_DFM (data=wvl, p=2, argvals=t, Lambdas=Lambdas, Basis=Basis,B = 500,alpha=0.05,r=0.8,cp=17:112)Del último conjunto de datos, conformado por las últimas 84 observaciones que presentan la misma media, se seleccionó el 80%, correspondiente a 76 observaciones, para la Fase I. Dentro de este conjunto, 37 observaciones se utilizaron para la estimación de los límites de control, mientras que las 20 observaciones restantes se destinaron a la Fase II, con el propósito de realizar el monitoreo del proceso.
A partir de la Fase I, el límite de control obtenido fue 2.163274
# Índices que NO se quieren considerar
indices_excluir <- c(
CL_DFM$indices,
CL_DFM$removed_indices
)
# Índices válidos
indices_validos <- setdiff(
1:nrow(Scores_g3),
indices_excluir
)
# Distancia funcional de Mahalanobis solo para índices válidos
DFM_final_3 <- sqrt(
rowSums(
(Scores_g3[indices_validos, ]^2) / matrix(
Eigenval_g3,
nrow = length(indices_validos),
ncol = length(Eigenval_g3),
byrow = TRUE
)
)
)
# Límite de control
CL_DFM_val <- CL_DFM$final_result$CL_DFM
# Data frame con índices reales
df_g3 <- data.frame(
obs = indices_validos,
DFM = DFM_final_3
)
limit_color <- "#D62728"
dfm_color <- "#1F77B4"
ggplot(df_g3, aes(x = obs, y = DFM)) +
geom_line(
color = dfm_color,
linewidth = 0.9
) +
geom_point(
color = dfm_color,
size = 1.6,
alpha = 0.85
) +
geom_hline(
yintercept = CL_DFM_val,
color = limit_color,
linetype = "dashed",
linewidth = 1
) +
scale_x_continuous(
labels = function(x) x + 16
) +
labs(
x = "Índice de observación",
y = "Valor de la DFM"
) +
theme_minimal(base_size = 13) +
theme(
panel.grid.minor = element_blank(),
plot.title = element_text(
face = "bold",
hjust = 0.5
)
)
| Observación real | Valor DFM |
|---|---|
| 104 | 2.202 |
| 109 | 2.598 |
Dado que el MFPCA se basa en la descomposición de Karhunen–Loève, se tiene que:
El vector funcional multivariante
\[\mathbf{X} = \big(X^{(1)}, \dots, X^{(p)}\big)\]
admite una representación finita de Karhunen–Loève si y solo si cada uno de los procesos univariantes \(X^{(j)}, \quad j = 1, \dots, p\) admite una representación finita de Karhunen–Loève.
Esto implica que cada proceso univariante pertenece a un espacio de Hilbert, en particular al espacio de funciones aleatorias cuadrado integrables \(L^2(\mathcal{T})\), y que el operador de covarianza asociado es de tipo Hilbert–Schmidt, autoadjunto (simétrico) y definido positivo.
Además, otro supuesto importante es que la función media multivariada sea constante, es decir, se quiere juzgar el sistema de hipótesis
\[ H_0: \mathbb{E}[\mathbf{X}_1(t)] = \cdots = \mathbb{E}[\mathbf{X}_N(t)] \] donde
\[ \mathbb{E}[\mathbf{X}_i(t)] = \boldsymbol{\mu}_i(t) = \left( \mu_i^{(1)}(t), \mu_i^{(2)}(t), \ldots, \mu_i^{(p)}(t) \right) \]
es la función media multivariada de la observación (i).
Bajo la hipótesis alternativa, se asume la existencia de al menos un punto de cambio (k^{}), tal que la media funcional multivariada cambia antes y después de dicho punto. Por tanto, el modelo puede escribirse como
\[ X_i(t) = \begin{cases} \mu^{(1)}(t) + \varepsilon_i(t), & 1 \leq i \leq k^{\ast}, \\[6pt] \mu^{(2)}(t) + \varepsilon_i(t), & k^{\ast} < i \leq N, \end{cases} \]
donde (^{(1)}(t)) y (^{(2)}(t)) son funciones medias multivariadas, y (_i(t)) representa un término de error.
Para reducir la dimensión del problema, se aplica un análisis de componentes principales funcionales multivariado. Primero se estima la media funcional multivariada muestral como
\[ \bar{\mathbf{X}}_N(t) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{X}_i(t). \]
Luego, cada observación centrada se aproxima mediante la expansión
\[ \mathbf{X}_i(t) - \bar{\mathbf{X}}_N(t) \approx \sum_{\ell=1}^{d} \eta_{i,\ell}\boldsymbol{\psi}_{\ell}(t), \]
donde \(d\) es el número de componentes principales retenidas, \(\boldsymbol{\psi}_{\ell}(t)\) es la \(\ell\)-ésima función propia multivariada y \(\eta_{i,\ell}\) es el score asociado a la observación \(i\) en la componente \(\ell\).
Cada función propia multivariada tiene la forma
\[ \boldsymbol{\psi}_{\ell}(t) = \left( \psi_{\ell}^{(1)}(t), \psi_{\ell}^{(2)}(t), \ldots, \psi_{\ell}^{(p)}(t) \right). \] En particular,
\[ \mathcal{H} = L^2(\mathcal{T}_1) \times L^2(\mathcal{T}_2) \times \cdots \times L^2(\mathcal{T}_p). \]
El producto interno en este espacio puede definirse como
\[ \left\langle f,g \right\rangle_{\mathcal{H}} = \sum_{j=1}^{p} \int_{\mathcal{T}_j} f^{(j)}(t)g^{(j)}(t)\,dt. \]
A partir de los scores de la MFPCA, se construye un estadístico acumulativo. Para cada posible punto de cambio \(k\), se define
\[ T_N(k) = \frac{1}{N} \sum_{\ell=1}^{d} \frac{1}{\hat{\lambda}_{\ell}} \left( \sum_{i=1}^{k} \hat{\eta}_{i,\ell} - \frac{k}{N} \sum_{i=1}^{N} \hat{\eta}_{i,\ell} \right)^2. \]
donde \(\hat{\lambda}_{\ell}\) es el valor propio asociado a la componente principal funcional multivariada \(\ell\), y \(\hat{\eta}_{i,\ell}\) son los scores estimados.
El estadístico global de prueba se obtiene acumulando la evidencia de cambio a lo largo de todos los posibles puntos \(k\):
\[ S_{N,d} = \frac{1}{N} \sum_{k=1}^{N} T_N\left(k\right). \]
Equivalentemente,
\[ S_{N,d} = \frac{1}{N^2} \sum_{\ell=1}^{d} \frac{1}{\hat{\lambda}_{\ell}} \sum_{k=1}^{N} \left( \sum_{i=1}^{k} \hat{\eta}_{i,\ell} - \frac{k}{N} \sum_{i=1}^{N} \hat{\eta}_{i,\ell} \right)^2. \]
El punto de cambio estimado se obtiene como el valor de \(k\) que maximiza el estadístico:
\[ \hat{\theta}_N = \inf \left\{ x: T_N(x) = \sup_{0 \leq y \leq 1} T_N(y) \right\}. \]
Para verificar los supuestos se crea el objeto con ayuda de la función multiFunData().
p <- length(fit)
X_list <- lapply(1:p, function(j) {
fd2funData(fdobj = fit[[j]]$fd, argvals = grid)
})
mFData <- multiFunData(X_list)Se usa la metodología mencionada anteriormente para verificar si la media es constante.
choose_d <- function(lambda, threshold = 0.90) {
pve <- lambda / sum(lambda)
d <- min(which(cumsum(pve) >= threshold))
return(d)
}
critical_Kd <- function(d, alpha = 0.05, B = 20000, J = 500, seed = 123) {
set.seed(seed)
r <- seq_len(J)
weights <- rep(1 / (pi^2 * r^2), times = d)
Z <- matrix(rnorm(B * d * J), nrow = B)
K_sim <- drop((Z^2) %*% weights)
unname(quantile(K_sim, probs = 1 - alpha))
}
get_mfpca_segment <- function(fits, fit, grid, idx,
M_max = 10,
K_uni = 10,
threshold = 0.90) {
p <- dim(fits)[3]
fd_list <- lapply(1:p, function(j) {
mat <- fits[, idx, j, drop = FALSE]
mat <- matrix(mat, nrow = length(grid), ncol = length(idx))
funData(
argvals = grid,
X = t(mat)
)
})
mFData <- multiFunData(fd_list)
pca_list <- lapply(1:p, function(j) {
fdobj <- fit[[j]]$fd
fd_segment <- fdobj
fd_segment$coefs <- fdobj$coefs[, idx, drop = FALSE]
pca.fd(
fdobj = fd_segment,
nharm = K_uni,
centerfns = TRUE
)
})
eigenfun.funData <- lapply(1:p, function(j) {
fd2funData(
pca_list[[j]]$harmonics,
argvals = grid
)
})
scores_list <- lapply(1:p, function(j) {
pca_list[[j]]$scores
})
uniExpansions <- lapply(1:p, function(j) {
list(
type = "given",
functions = eigenfun.funData[[j]],
scores = scores_list[[j]],
ortho = TRUE
)
})
MFPCAObj <- MFPCA(
mFData = mFData,
M = M_max,
uniExpansions = uniExpansions,
fit = FALSE,
verbose = FALSE
)
lambda <- MFPCAObj$values
d <- choose_d(lambda, threshold = threshold)
list(
scores = MFPCAObj$scores[, 1:d, drop = FALSE],
lambda = lambda[1:d],
d = d,
pve = cumsum(lambda / sum(lambda)),
MFPCAObj = MFPCAObj
)
}
test_segment_local_mfpca <- function(fits, fit, grid, start, end,
M_max = 10,
K_uni = 10,
threshold = 0.90,
alpha = 0.05,
min_size = 5) {
idx <- start:end
n <- length(idx)
if (n < 2 * min_size) {
return(list(
reject = FALSE,
statistic = NA,
critical_value = NA,
d = NA,
k_hat = NA,
k_local = NA,
T_k = rep(NA, n)
))
}
mf <- get_mfpca_segment(
fits = fits,
fit = fit,
grid = grid,
idx = idx,
M_max = M_max,
K_uni = K_uni,
threshold = threshold
)
scores <- as.matrix(mf$scores)
lambda <- mf$lambda
d <- mf$d
C <- matrix(0, nrow = n, ncol = d)
for (m in 1:d) {
C[, m] <- cumsum(scores[, m]) -
(seq_len(n) / n) * sum(scores[, m])
}
T_k <- rowSums(sweep(C^2, 2, lambda, "/")) / n
S <- sum(T_k) / n
critical_value <- critical_Kd(
d = d,
alpha = alpha
)
valid_k <- min_size:(n - min_size)
k_local <- valid_k[which.max(T_k[valid_k])]
k_global <- start + k_local - 1
list(
reject = S > critical_value,
statistic = S,
critical_value = critical_value,
d = d,
k_hat = k_global,
k_local = k_local,
T_k = T_k,
pve = mf$pve
)
}
MfChange <- function(fits, fit, grid,
M_max = 10, K_uni = 10,
threshold = 0.90, alpha = 0.05,
min_size = 30) {
N <- dim(fits)[2]
change_points <- c()
results <- list()
search_segment <- function(start, end) {
test <- test_segment_local_mfpca(
fits = fits,
fit = fit,
grid = grid,
start = start,
end = end,
M_max = M_max,
K_uni = K_uni,
threshold = threshold,
alpha = alpha,
min_size = min_size
)
results[[length(results) + 1]] <<- data.frame(
start = start,
end = end,
n = end - start + 1,
d = test$d,
statistic = test$statistic,
critical_value = test$critical_value,
reject = test$reject,
k_hat = test$k_hat
)
if (!test$reject || is.na(test$k_hat)) {
return(NULL)
}
k <- test$k_hat
change_points <<- c(change_points, k)
search_segment(start, k)
search_segment(k + 1, end)
}
search_segment(1, N)
list(
change_points = sort(change_points),
results = do.call(rbind, results)
)
}ChangeDetectionMulti <- MfChange(
fits = fits,
fit = fit,
grid = grid,
M_max = 10,
K_uni = 10,
threshold = 0.90,
alpha = 0.05,
min_size = 5
)Hay evidencia estadísticamente significativa para rechazar $H_0: [_1(t)] = = [_N(t)]. Se muestran los resultados a continuación en la Tabla 3:
load("ChangeDetectionMulti.RData")
separar_cp_mfd <- function(mFData, puntos_cambio, na.rm = FALSE) {
N <- funData::nObs(mFData)
puntos <- sort(unique(as.integer(puntos_cambio)))
inicios <- c(1, puntos + 1)
fines <- c(puntos, N)
lista_resultados <- list()
for (i in seq_along(inicios)) {
idx <- inicios[i]:fines[i]
mfd_sub <- funData::extractObs(mFData, obs = idx)
mfd_centrado <- mfd_sub - funData::meanFunction(mfd_sub, na.rm = na.rm)
nombre_grupo <- paste0("Segmento_", i)
lista_resultados[[nombre_grupo]] <- list(
n = length(idx),
indices = idx,
rango = paste(inicios[i], "a", fines[i]),
mFData_original = mfd_sub,
mFData_centrado = mfd_centrado,
media = funData::meanFunction(mfd_sub, na.rm = na.rm)
)
}
return(lista_resultados)
}
cpMulti <- separar_cp_mfd(mFData, ChangeDetectionMulti$change_points)
ns <- unlist(lapply(1:length(cpMulti), function(x){cpMulti[[x]]$n}))
rangoss <- unlist(lapply(1:length(cpMulti), function(x){cpMulti[[x]]$rango}))
cptable <- data.frame(Segmento = 1:length(cpMulti),
n = ns,
Rango = rangoss)
cptable %>%
kable(format = "pandoc",
digits = 4,
align = "c",
escape = FALSE)| Segmento | n | Rango |
|---|---|---|
| 1 | 84 | 1 a 84 |
| 2 | 12 | 85 a 96 |
| 3 | 15 | 97 a 111 |
| 4 | 23 | 112 a 134 |
| 5 | 5 | 135 a 139 |
| 6 | 5 | 140 a 144 |
| 7 | 11 | 145 a 155 |
| 8 | 15 | 156 a 170 |
| 9 | 5 | 171 a 175 |
| 10 | 5 | 176 a 180 |
| 11 | 16 | 181 a 196 |
| 12 | 7 | 197 a 203 |
| 13 | 39 | 204 a 242 |
| 14 | 14 | 243 a 256 |
| 15 | 7 | 257 a 263 |
| 16 | 5 | 264 a 268 |
En las siguientes figuras se muestran las curvas y funciones media para el segmento 1 y segmento 13.
Para determinar las componentes principales funcionales multivariadas, se seguirá el enfoque propuesto por (Happ y Greven 2018), el cual permite realizar MFPCA incluso cuando los dominios de los procesos funcionales son diferentes.
En el presente caso, los dominios son iguales; no obstante, este enfoque resulta conveniente, ya que garantiza propiedades deseables del operador de covarianza del proceso funcional multivariado, tales como ser de tipo Hilbert–Schmidt, autoadjunto y definido positivo.
Dadas muestras \(\mathbf{\chi_1}, \dots, \mathbf{\chi_N}\) de \(\mathbf{\chi}\), el procedimiento de estimación propuesto para MFPCA consta de cuatro pasos:
X.G1_list <- lapply(1:p, function(j){
funData2fd(cpMulti[[idxSegMulti[1]]]$mFData_original[[j]],
basisobj = Basis,
lambda = Lambdas[j])
})
Mj <- 2
pca.G1_list <- lapply(1:p, function(j){
pca.fd(fdobj = X.G1_list[[j]],
nharm = Mj,
centerfns = TRUE)
})| Proceso | Varianza_Explicada |
|---|---|
| \(\chi^{(1)}\) | 0.9928 |
| \(\chi^{(2)}\) | 0.9905 |
| \(\chi^{(3)}\) | 0.9925 |
| \(\chi^{(4)}\) | 0.9706 |
| \(\chi^{(5)}\) | 0.9948 |
| \(\chi^{(6)}\) | 0.9961 |
| \(\chi^{(7)}\) | 0.9929 |
Definir la matriz \(\Xi \in \mathbb{R}^{N \times M_+}\), donde cada fila \(\Xi_i\) contiene todos los scores estimados para una observación: \(\Xi_i = \big(\hat{\xi}^{(1)}_{i,1}, \dots, \hat{\xi}^{(1)}_{i,M_1}, \dots, \hat{\xi}^{(p)}_{i,1}, \dots, \hat{\xi}^{(p)}_{i,M_p}\big).\) Un estimador \(\hat{Z} \in \mathbb{R}^{M_+ \times M_+}\) de la matriz de bloques \(Z\) se obtiene como \(\hat{Z} = (N - 1)^{-1} \Xi^\top \Xi.\)
Realizar un análisis espectral de matrices para \(\hat{Z}\), obteniendo valores propios \(\hat{\nu}_m\) y vectores propios ortonormales \(\hat{\mathbf{c}}_m\), \(m = 1, \dots, M_+\).
Los estimadores elemento a elemento de las funciones propias multivariantes asociadas a \(\hat{\nu}_m\) están dados por: \(\hat{\psi}^{(j)}_m(t_j) = \sum_{n=1}^{M_j} [\hat{\mathbf{c}}_m]^{(j)}_n \, \hat{\phi}^{(j)}_n(t_j), t_j \in \mathcal{T}_j,\) y los scores multivariantes pueden calcularse como:
\[\rho_m = \sum_{j=1}^{p} \sum_{n=1}^{M_j} [\boldsymbol{c_m}]_{n}^{(j)} \xi_{n}^{(j)}\] Dado que se identificaron múltiples puntos de cambio, para el desarrollo de los siguientes análisis se seleccionará el segmento 1, el cual contiene 84 observaciones funcionales multivariadas
eigenfun.G1_fD <- lapply(1:p, function(j){
fd2funData(pca.G1_list[[j]]$harmonics, argvals = grid)
})
scores.G1_list <- lapply(1:p, function(j){
pca.G1_list[[j]]$scores
})
uniExpansions <- lapply(1:p, function(j){
list(type = "given",
functions = eigenfun.G1_fD[[j]],
scores = scores.G1_list[[j]])
})
MFPCA.G1 <- MFPCA(cpMulti[[idxSegMulti[1]]]$mFData_original,
M = 3,
uniExpansions = uniExpansions,
fit = TRUE,
verbose = FALSE)
\[\text{Cov}(\mathbf{X}(t),\mathbf{X}(s))\]
\[\begin{pmatrix} \Sigma_{11} & \Sigma_{12} & \Sigma_{13} & \Sigma_{14} & \Sigma_{15} & \Sigma_{16} & \Sigma_{17} \\ \Sigma_{21} & \Sigma_{22} & \Sigma_{23} & \Sigma_{24} & \Sigma_{25} & \Sigma_{26} & \Sigma_{27} \\ \Sigma_{31} & \Sigma_{32} & \Sigma_{33} & \Sigma_{34} & \Sigma_{35} & \Sigma_{36} & \Sigma_{37} \\ \Sigma_{41} & \Sigma_{42} & \Sigma_{43} & \Sigma_{44} & \Sigma_{45} & \Sigma_{46} & \Sigma_{47} \\ \Sigma_{51} & \Sigma_{52} & \Sigma_{53} & \Sigma_{54} & \Sigma_{55} & \Sigma_{56} & \Sigma_{57} \\ \Sigma_{61} & \Sigma_{62} & \Sigma_{63} & \Sigma_{64} & \Sigma_{65} & \Sigma_{66} & \Sigma_{67} \\ \Sigma_{71} & \Sigma_{72} & \Sigma_{73} & \Sigma_{74} & \Sigma_{75} & \Sigma_{76} & \Sigma_{77} \end{pmatrix}\]
\(\Sigma_{\chi^{(1)}\chi^{(1)}}(t, s) = (N - 1)^{-1} \sum_{i=1}^{N} \{x_{i}^{(1)}(t) - \bar{x}^{(1)}(t)\}\{x_{i}^{(1)}(s) - \bar{x}^{(1)}(s)\}\)
\(\Sigma_{\chi^{(1)}\chi^{(2)}}(t, s) = (N - 1)^{-1} \sum_{i=1}^{N} \{x_{i}^{(1)}(t) - \bar{x}^{(1)}(t)\}\{x_{i}^{(2)}(s) - \bar{x}^{(2)}(s)\}\)