Taller 2 - Análisis de Datos Funcionales

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

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

1 1. FPCA: Aplicación

1.1 a. Seleccione uno de los 7 procesos que hacen parte del estudio espectrométrico, y con las observaciones muestrales del proceso seleccionado:

1.1.1 1) Supuestos del FPCA

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 el proceso se encuentra centrado. Por lo tanto, a cada observación se le restará la media estimada, como se muestra a continuación en la Figura 1.

X <- t(fits[, , 2])

X.fd <- fit[[2]]

grid <- seq(min(t), max(t), length = 1000)

mean <- colMeans(X)

Xc <- sweep(X, 2, mean, FUN = "-")

mean.Xc <- colMeans(Xc)

Figura 1: Curvas espectrométricas centradas junto con su media

1.1.2 2) Funciones Propias

Luego de haber centrado el proceso, se encuentran las estimaciones de las funciones propias asociadas al operador de covarianza del proceso funcional.

pcaObj <- pca.fd(X.fd$fd, nharm = 3, centerfns = TRUE)

eigenfn <- eval.fd(grid, pcaObj$harmonics)

eigenval <- pcaObj$values[1:3]

Las tres primeras funciones propias se pueden visualizar en la Figura 2.

[1] "done"

Figura 2: Gráfico de las tres primeras funciones propias estimadas

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 3, Figura 4 y Figura 5.

Figura 3: Efecto de las primera componente principal sobre la media del proceso

Figura 4: Efecto de lsegunda componente principal sobre la media del proceso

Figura 5: Efecto de la tercera componente principal sobre la media del proceso

• 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.

1.1.3 3) Varianza explicada

En las gráficas anteriores es posible ver cuánto porcentaje de variabilidad es explicado por las tres funciones propias estimadas. Esto también se muestra a través de un scree plot en la Figura 6.

Figura 6: Scree plot de las tres primeras componentes principales

De esta forma, es posible ver que las dos primeras componentes principales ya explican más del 90% de variabilidad total.

1.1.4 4) Descomposición Karhunen-Loève

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 el 96.33 de la variabilidad, aproximadamente. Por lo tanto, para cada observación funcional, se tiene su aproximación \[ X_i(t) \approx \hat \mu(t) + \sum_{k=1}^K \hat {\xi}_{ik} \hat \nu _ k(t) \] con \(i = 1,\ldots,268\) y \(K=2\).

A continuación, se expresan las tres primeras observaciones muestrales con la descomposición de Karhunen-Loève. En la @ se pueden ver las observaciones junto con su aproximación.

Scores <- matrix(pcaObj$scores[,1:2], nrow = 268, ncol = 2)

eigenfns <- eval.fd(evalarg = grid, fdobj = pcaObj$harmonics)[,1:2]

mu <- eval.fd(evalarg = grid, fdobj = pcaObj$meanfd)

X.approx <- Scores %*% t(eigenfns)
X.approx <- sweep(X.approx, 2, mu, "+") 

Figura 7: Curvas observadas y sus aproximaciones mediante la descomposición de Karhunen-Loève con las dos primeras funciones propias

Figura 8: Curvas observadas y sus aproximaciones mediante la descomposición de Karhunen-Loève con las dos primeras funciones propias

Figura 9: Curvas observadas y sus aproximaciones mediante la descomposición de Karhunen-Loève con las dos primeras funciones propias

A continuación, en la Tabla 1 se muestra el error cuadrático medio de predicción para cada observación.

Tabla 1: Error cuadrático medio de las aproximaciones con la descomposición de Karhunen-Loève

	MPSE
Primera Observación	1.2944
Segunda Observación	3.8829
Tercera Observación	7.1369

1.1.5 5) Realizaciones del primer score

A continuación se presentan las realizaciones del primer score en la Figura 10.

Figura 10: Gráficos de las realizaciones del primer score

Es posible ver que los scores tienen una distribución asimétrica y también presenta algunos outliers por valores muy altos.

1.1.6 6) Gráficos de Dispersión de los Scores

A continuación, se muestra el gráfico de dispersión de las dos primeras componentes principales, que explican el 96.33 de la variabilidad, aproximadamente.

Figura 11: Score plot de las dos primeras componentes principales

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.

1.1.7 7) Función de Covarianza y Funciones Propias

En la Figura 12 se presenta la estimación de la función de covarianza y las funciones propias.

cov.fd <- var.fd(X.fd$fd)
cov.grid <- seq(min(t), max(t), length = 20)
cov.mat <- eval.bifd(cov.grid, cov.grid, cov.fd)

Figura 12: Gráfico de la función de covarianza estimada y las funciones propias estimadas

Como es posible visualizar, la forma de la primera función propia es similar a las que se pueden ver en la función de covarianza debido a su alto porcentaje de variabilidad explicada.

1.1.8 8) Teorema de Mercer

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

En la Figura 13 se puede ver la función de covarianza junto con su aproximación y en la Figura 14 se ven sus residuales.

eigenfn2 <- eval.fd(cov.grid, pcaObj$harmonics)

NU <- as.matrix(eigenfn2)

D <- diag(eigenval)

mercer <- NU %*% D %*% t(NU)

cov.diff <- mercer - cov.mat

Figura 13: Gráfico de la función de covarianza y su aproximación por el teorema de Mercer

Figura 14: Gráfico de la diferencia entre la función de covarianza y su aproximación por el teorema de Mercer

1.2 b. Cartas de Control

1.2.1 FPCA

1.2.2 Distancia Funcional de Mahalanobis

2 2. MFPCA

2.1 Considere todo los 7 procesos espectrométricos del proceso de producción de azúcar:

2.1.1 Cuáles son los supuestos que requiere el MFPCA

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.

2.1.2 Verique que estos supuestos se cumplen

3 Usando MFPCA:

3.1 Encuentre las estimaciones de las tres primeras funciones propias multivariadas, grafíquelas y presente una explicación de ellas.

Para determinar las componentes principales funcionales multivariadas, se seguirá el enfoque propuesto por Clara Happ, 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.

Xc1 <- center.fd(fit[[1]]$fd)
Xc2 <- center.fd(fit[[2]]$fd)
Xc3 <- center.fd(fit[[3]]$fd)
Xc4 <- center.fd(fit[[4]]$fd)
Xc5 <- center.fd(fit[[5]]$fd)
Xc6 <- center.fd(fit[[6]]$fd)
Xc7 <- center.fd(fit[[7]]$fd)

Dadas muestras \(\mathbf{\chi_1}, \dots, \mathbf{\chi_N}\) de \(\mathbf{\chi}\), el procedimiento de estimación propuesto para MFPCA consta de cuatro pasos:

1. Para cada elemento \(\chi^{(j)}\) se estima un FPCA univariante basado en las observaciones \(\chi^{(j)}_1, \dots, \chi^{(j)}_N\). Esto produce funciones propias estimadas \(\hat{\phi}^{(j)}_m\) y scores \(\hat{\xi}^{(j)}_{i,m}\), \(i = 1, \dots, N\), \(m = 1, \dots, M_j\),para soluciones de truncamiento elegidas adecuadamente \(M_j\).

pcaObj1 <- pca.fd(Xc1, nharm = 2)
pcaObj2 <- pca.fd(Xc2, nharm = 2)
pcaObj3 <- pca.fd(Xc3, nharm = 2)
pcaObj4 <- pca.fd(Xc4, nharm = 2)
pcaObj5 <- pca.fd(Xc5, nharm = 2)
pcaObj6 <- pca.fd(Xc6, nharm = 2)
pcaObj7 <- pca.fd(Xc7, nharm = 2)

pca_list <- list(pcaObj1, pcaObj2, pcaObj3, pcaObj4, 
                 pcaObj5, pcaObj6, pcaObj7)

var_exp <- sapply(pca_list, function(obj) {
  sum(obj$varprop)
})

tabla <- data.frame(
  Proceso = paste0("$\\chi^{(", 1:7, ")}$"),
  Varianza_Explicada = round(var_exp, 4)
)

knitr::kable(tabla, escape = FALSE)

Proceso	Varianza_Explicada
\(\chi^{(1)}\)	0.9424
\(\chi^{(2)}\)	0.9633
\(\chi^{(3)}\)	0.9580
\(\chi^{(4)}\)	0.9645
\(\chi^{(5)}\)	0.9539
\(\chi^{(6)}\)	0.9416
\(\chi^{(7)}\)	0.9483

2. 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.\)

3. 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_+\).

4. 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:

Xc_list <- list(Xc1, Xc2, Xc3, Xc4, Xc5, Xc6, Xc7)
fit.funData <- list()

for (i in 1:7) {
  fit.funData[[i]] <- fd2funData(Xc_list[[i]], argvals = grid)
}

mFData <- multiFunData(fit.funData)

uniExpansions <- list(list(type = "uFPCA", npc = 2), 
                      list(type = "uFPCA", npc = 2),
                      list(type = "uFPCA", npc = 2), 
                      list(type = "uFPCA", npc = 2),
                      list(type = "uFPCA", npc = 2),
                      list(type = "uFPCA", npc = 2),
                      list(type = "uFPCA", npc = 2))

MFPCAObj <- MFPCA(mFData, 
                  M = 3, 
                  uniExpansions = uniExpansions, 
                  fit = TRUE,
                  verbose = TRUE)
save(MFPCAObj, file = "MFPCAObj.RData")

load("MFPCAObj.RData")

par(mfrow = c(2,4), mar = c(3,3,2,1))  

cols <- 1:3

for (j in 1:7) {
  
  matplot(MFPCAObj$functions[[j]]@argvals[[1]],
          t(MFPCAObj$functions[[j]]@X),
          type = "l", lty = 1,
          col = cols,
          xlab = "t", ylab = "",
          main = paste("Proceso", j))
  
  legend("topright",
         legend = c(expression(psi[1]),
                    expression(psi[2]),
                    expression(psi[3])),
         col = cols, lty = 1, lwd = 2,
         cex = 0.8, bty = "n")
}

Xc_list <- list(Xc1, Xc2, Xc3, Xc4, Xc5, Xc6, Xc7)
par(mfrow = c(2,4), mar = c(3,3,2,1))  # 7 plots (2x4)
for (j in 1:7) {
  plot(Xc_list[[j]],
          xlab = "t", ylab = "",
          main = paste("Proceso", j))
}

3.2 Grafique Σ11 y Σ12.

\(\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)\}\)

wvl_c1 <- sweep(wvl[,,1], 1, rowMeans(wvl[,,1]), "-")
wvl_c2 <- sweep(wvl[,,2], 1, rowMeans(wvl[,,2]), "-")

sigma11 = (wvl_c1 %*% t(wvl_c1))/(268-1)
sigma12 = (wvl_c1 %*% t(wvl_c2))/(268-1)

escala_colores <- list(
  list(0, "#000033"),      
  list(0.02, "#0000FF"),   
  list(0.05, "#FFFFFF"),   
  list(0.1, "#FFFF00"),    
  list(0.3, "#FF7F00"),    
  list(0.6, "#FF0000"),    
  list(1, "#4B0082")       
)

p <- plot_ly(x = ~t, y = ~t, z = ~sigma11) %>% 
  add_surface(
    colorscale = escala_colores,
    colorbar = list(title = "Covarianza/"),
    contours = list(
      z = list(
        show = TRUE, 
        usecolormap = TRUE, 
        highlightcolor = "#ff0000", 
        project = list(z = TRUE) 
      )
    )
  ) %>%
  layout(
    title = list(text = "Superficie de Covarianza Σ11"),
    scene = list(
      xaxis = list(title = "Longitud de onda (nm)"),
      yaxis = list(title = "Longitud de onda (nm)"),
      zaxis = list(
        title = "Valor", 
        range = c(0,4567),
        tickformat = ".0f"
      )
    )
  )
p

\(\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)\}\)

escala_colores <- list(
  list(0, "#000033"),      
  list(0.02, "#0000FF"),   
  list(0.05, "#FFFFFF"),   
  list(0.1, "#FFFF00"),    
  list(0.3, "#FF7F00"),    
  list(0.6, "#FF0000"),    
  list(1, "#4B0082")       
)

p <- plot_ly(x = ~t, y = ~t, z = ~sigma11) %>% 
  add_surface(
    colorscale = escala_colores,
    colorbar = list(title = "Covarianza/"),
    contours = list(
      z = list(
        show = TRUE, 
        usecolormap = TRUE, 
        highlightcolor = "#ff0000", 
        project = list(z = TRUE) 
      )
    )
  ) %>%
  layout(
    title = list(text = "Superficie de Covarianza Σ12"),
    scene = list(
      xaxis = list(title = "Longitud de onda (nm)"),
      yaxis = list(title = "Longitud de onda (nm)"),
      zaxis = list(
        title = "Valor", 
        range = c(0,827),
        tickformat = ".0f"
      )
    )
  )
p

var_fd <- var.fd(fit[[1]]$fd)
Cmat <- eval.bifd(t, t, var_fd)
plot(t,diag(Cmat))

3.3 Cuántas funciones propias del operador explican el 90 % de la varianza total del proceso. Grafíquelas.

var_exp <- MFPCAObj14[["values"]] / sum(MFPCAObj14[["values"]])

cum_var_exp <- cumsum(var_exp)

plot(cum_var_exp[1:8], type = "b",
     xlab = "Número de componentes",
     ylab = "Varianza explicada acumulada",
     ylim = c(0, 1.2),
     pch = 19)

abline(h = 0.9, col = "red", lty = 2)

text(1:length(cum_var_exp),
     cum_var_exp,
     labels = round(cum_var_exp, 3),
     pos = 3,   
     cex = 0.8)

load("MFPCAObj2.RData")

par(mfrow = c(2,4), mar = c(3,3,2,1))  

cols <- 1:2

for (j in 1:7) {
  
  matplot(MFPCAObj2$functions[[j]]@argvals[[1]],
          t(MFPCAObj2$functions[[j]]@X),
          type = "l", lty = 1,
          col = cols,
          xlab = "t", ylab = "",
          main = paste("Proceso", j))
  
  legend("topright",
         legend = c(expression(psi[1]),
                    expression(psi[2])),
         col = cols, lty = 1, lwd = 2,
         cex = 0.8, bty = "n")
}

3.4 Grafique los valores propios del MFPCA y explique su comportamiento.

values <- MFPCAObj14[["values"]]

plot(log(values), type = "b",
     pch = 19, lwd = 2,
     col = "darkred",
     ylim = c(0,14),
     xlab = "Número de componentes",
     ylab = "Valor propio (escala log)")

grid()

3.5 Presente una tabla con las correlaciones de los scores asociados a los dos primeros procesos.

round(colMeans(pcaObj1$scores),4)

[1] 0 0

round(colMeans(pcaObj2$scores),4)

[1] 0 0

cor_matrix <- function(obj1, obj2){
  matrix(c(
    cor(obj1$scores[,1], obj2$scores[,1]),
    cor(obj1$scores[,1], obj2$scores[,2]),
    cor(obj1$scores[,2], obj2$scores[,1]),
    cor(obj1$scores[,2], obj2$scores[,2])
  ), nrow = 2, byrow = TRUE)
}


cor_scores_12 <- cor_matrix(pcaObj1,pcaObj2)

rownames(cor_scores_12) <- c("E11", "E21")
colnames(cor_scores_12) <- c("E12", "E22")

knitr::kable(cor_scores_12, digits = 3)

	E12	E22
E11	0.800	0.405
E21	0.566	-0.329

3.6 Grafique las correlación entre los scores multivariados usando un mapa de calor.

pca_list <- list(pcaObj1, pcaObj2, pcaObj3, pcaObj4, 
                 pcaObj5, pcaObj6, pcaObj7)

n <- length(pca_list)

cor_list <- list()

for(i in 1:n){
  for(j in i:n){
    cor_list[[paste0("Proc",i,"-Proc",j)]] <- cor_matrix(pca_list[[i]], pca_list[[j]])
  }
}

# inicializar matriz grande (2 PCs por proceso)
cor_big <- matrix(NA, nrow = 2*n, ncol = 2*n)

# nombres
labels <- unlist(lapply(1:n, function(i) paste0("P", i, "_PC", 1:2)))
rownames(cor_big) <- labels
colnames(cor_big) <- labels

for(i in 1:n){
  for(j in 1:n){
    
    # matriz 2x2 entre proceso i y j
    block <- cor_matrix(pca_list[[i]], pca_list[[j]])
    
    # posiciones en la matriz grande
    rows <- ((i-1)*2 + 1):(i*2)
    cols <- ((j-1)*2 + 1):(j*2)
    
    cor_big[rows, cols] <- round(block,2)
  }
}

cor_long <- melt(cor_big)

ggplot(cor_long, aes(Var1, Var2, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white",
                       midpoint = 0, limits = c(-1,1)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "", fill = "Correlación")

cor_scores_multi <- cor(MFPCAObj14$scores)
cor_long <- melt(cor_scores_multi)

ggplot(cor_long, aes(Var1, Var2, fill = value)) +
  geom_tile() +
  geom_text(aes(label = round(value, 2)), size = 3) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white",
                       midpoint = 0, limits = c(-1,1)) +
  theme_minimal() +
  labs(x = "Componentes", y = "Componentes", fill = "Correlación") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

3.7 Presente una matriz con los graficos de dispersión de las realizaciones de los scores que conjuntamente explican el 90 % de la variación del proceso. Explique el comportamiento de los scores

df <- data.frame(
  PC1 = MFPCAObj2$scores[,1],
  PC2 = MFPCAObj2$scores[,2]
)

ggplot(df, aes(x = PC1, y = PC2)) +
  geom_point(alpha = 0.6, size = 2, color = "steelblue") +
  theme_minimal() +
  labs(
    title = "Dispersión de los scores (PC1 vs PC2)",
    x = "Componente principal 1",
    y = "Componente principal 2"
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray")

4 Exprese las tres primeras observaciones muestrales usando la descomposición de de Karhunen-Loève

En el caso multivariado, se tiene la descomposición de Karhunen–Loève, de forma que para el vector de curvas aleatorias, $ L^2( _1 ) L^2( _p) $, se tiene que

\[ \mathbf{X}(\mathbf{t}) = \boldsymbol{\mu}(\mathbf{t}) + \sum_{k=1}^\infty \xi_k \boldsymbol{\psi}_k (\mathbf{t}) \] con \(\mathbf{t} = (t_1,\ldots,t_p)^t\). En este caso, se tiene un vector de funciones de tamaño \(p=7\), donde el dominio es el mismo para cada proceso, \(\mathcal{T_1}=\cdots=\mathcal{T_7} = [275, 560]\). De esta forma, se aproxima cada una de las observaciones muestrales de la siguiente manera:

\[ \mathbf{X}_i(\mathbf{t}) \approx \hat {\boldsymbol{\mu}}(\mathbf{t}) + \sum_{k=1}^K \hat{\xi}_{ik} \hat{\boldsymbol{\psi}}_k (\mathbf{t}) \]

4.1 Grafique las tres primeras observaciones muestrales Vs las aproximaciones obtenidas por la descomposición de Karhunen-Loève

X_list <- list(Xc1, Xc2, Xc3, Xc4, Xc5, Xc6, Xc7)
n_proc <- length(X_list)

for(j in 1:n_proc){
  
  # aproximaciones
  aprox3 <- MFPCAObj2[["fit"]][[j]]@X[1:3,]
  
  # dominio
  t <- MFPCAObj2[["functions"]][[j]]@argvals[[1]]
  
  # funciones originales
  muestra3 <- t(eval.fd(t, X_list[[j]])[,1:3])
  
  # gráfico
  par(mfrow = c(1,3))
  
  for(i in 1:3){
    plot(t, muestra3[i,], type = "l", lwd = 2,
         col = "blue",
         ylim = range(c(muestra3, aprox3)),
         xlab = "t", ylab = "Valor",
         main = paste("Proceso", j, "- Curva", i))
    
    lines(t, aprox3[i,], col = "red", lwd = 2, lty = 2)
    
    legend("topright",
           legend = c("Original", "Aprox"),
           col = c("blue", "red"),
           lty = c(1,2),
           bty = "n")
  }
}