#| 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)
}
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)
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"
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.
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.
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.
De esta forma, es posible ver que 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 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, "+") 
A continuación, en la Tabla 1 se muestra el error cuadrático medio de predicción para cada observación.
| MPSE | |
|---|---|
| Primera Observación | 1.2944 |
| Segunda Observación | 3.8829 |
| Tercera Observación | 7.1369 |
A continuación se presentan las realizaciones del primer score en la Figura 10.
Es posible ver que los scores tienen una distribución asimétrica y también presenta algunos outliers por valores muy altos.
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.
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.
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)
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.
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
T2.CL <- function(data, p, argvals, GCV.lambda, Basis, ncomp, alpha, B, seed){
X <- data[ , , p]
N <- dim(X)[2]
q <- numeric(B)
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = GCV.lambda[p])
set.seed(seed)
for (b in seq_len(B)) {
idx <- sample(1:N, size = N, replace = TRUE)
samples <- X[ , idx, drop = FALSE]
smoothlist <- smooth.basis(argvals, samples, fdParobj)
fdObj <- smoothlist$fd
pcaObj <- pca.fd(fdobj = fdObj, nharm = ncomp)
Scores <- as.matrix(pcaObj$scores)
Eigenval <- pcaObj$values[1:ncomp]
T2 <- diag(Scores %*% diag(1/Eigenval) %*% t(Scores))
q[b] <- max(T2)
}
CL.Boot <- quantile(q, probs = 1 - alpha)
smoothlist <- smooth.basis(argvals, X, fdParobj)
fdObj <- smoothlist$fd
pcaObj <- pca.fd(fdobj = fdObj, nharm = ncomp)
Scores <- as.matrix(pcaObj$scores)
Eigenval <- pcaObj$values[1:ncomp]
T2 <- diag(Scores %*% diag(1/Eigenval) %*% t(Scores))
return(list(MDens = q, CL = as.numeric(CL.Boot), T2 = T2))
}Muestras training y testing
ntotal <- nrow(X)
ntrain <- round(0.8 * ntotal)
ntest <- ntotal - ntrain
set.seed(123)
idx <- sample(1:ntotal, size = ntrain, replace = FALSE)
training <- wvl[ , idx, ]
testing <- wvl[ , -idx, ]Calculando el CL
CLobj <- T2.CL(data = training,
p = 2,
argvals = t,
GCV.lambda = Lambdas,
Basis = Basis,
ncomp = 2,
alpha = 0.05,
B = 1000,
seed = 123)T2.training <- CLobj$T2
CL <- CLobj$CL
OC <- T2.training >= CL
OCw <- which(OC)
col <- ifelse(OC == TRUE, "red", "black")
plot(1:ntrain, T2.training,
type = "o", col = col,
ylim = c(0, CL + 0.1))
abline(h = CL, lty = 2)
T2.test <- function(data, p, argvals, Basis, GCV.lambda, ncomp){
X <- data[ , , p]
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = GCV.lambda[p])
smoothlist <- smooth.basis(argvals, X, fdParobj)
fdObj <- smoothlist$fd
pcaObj <- pca.fd(fdobj = fdObj, nharm = ncomp)
Scores <- as.matrix(pcaObj$scores)
Eigenval <- pcaObj$values[1:ncomp]
T2 <- diag(Scores %*% diag(1/Eigenval) %*% t(Scores))
return(T2)
}T2.testing <- T2.test(data = testing,
p = 2,
argvals = t,
Basis = Basis,
GCV.lambda = Lambdas,
ncomp = 2)
OC.test <- T2.testing >= CL
OCw <- which(OC)
col <- ifelse(OC == TRUE, "red", "black")
plot(1:ntest, T2.testing,
type = "o", col = col,
ylim = c(0, CL + 0.1))
abline(h = CL, lty = 2)
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.
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:
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 |
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:
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)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))
}
\(\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"
)
)
)
pvar_fd <- var.fd(fit[[1]]$fd)
Cmat <- eval.bifd(t, t, var_fd)
plot(t,diag(Cmat))
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")
}
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()
round(colMeans(pcaObj1$scores),4)[1] 0 0round(colMeans(pcaObj2$scores),4)[1] 0 0cor_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 |
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) <- labelsfor(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))
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")
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}) \]
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")
}
}
En la sección 4 del artículo de Happ, C., & Greven, S. se presentan varias simulaciones para medir el desempeño de su propuesta. Se tienen tres escenarios:
En cada caso, se tienen 100 datasets con \(N=250\) observaciones: \[ \mathbf{x}_i(\mathbf{t}) = \sum_{m=1}^{M} \rho_{im}\psi_m(\mathbf{t}) + \boldsymbol{\varepsilon}_i(\mathbf{t}), \qquad \boldsymbol{\varepsilon}_i(\mathbf{t}) \overset{\mathrm{iid}}{\sim} \mathrm{N}_p(0, \sigma^2 \mathbf{I}) \] con \(\mathbf{t} \in \mathcal{T}\) e \(i=1,\ldots, N = 250\).
Además se tienen otras configuraciones:
\[ \nu_m^{\text{exp}} = \exp{(-(m+1)/2)}, \qquad \nu_m^{\text{lin}} = (M+1-m)/M \]
La precisión de las estimaciones se mide con el error relativo.
En la primera configuración se tienen en cuenta las primeras \(M=8\) funciones de la base de Fourier en el intervalo \([0,1]\) para ambos procesos funcionales que forman el proceso bivariado de forma que \(\mathcal{T_1}=\mathcal{T_2}=[0,1]\).
Las obsevaciones se muestrean de un grid equidistante de \(S_1=S_2=100\) puntos y se usan \(M_1=M_2=8\) componentes principales funcionales univariadas para estimar las \(M=8\) componentes principales funcionales multivariadas. El método propuesto por los autores se compara con el método implementado en fda, que se denotará como \(\mathrm{MFPCA_{RS}}\), donde se realiza un suavizamiento previo con \(K=15\) splines cúbicos. La función pca.fd(), a diferencia de MFPCA(), devuelve la estimación de dos scores por observación bivariada y función propia estimada, de forma que se suman para poder compararlos con la nueva metodología.
La simulación se realiza mediante la función simMultiFunData() del paquete funData.
Sim_Setting1 <- function(eValType = c("exponential", "linear"),
s2 = 0,
seed = NULL,
N = 250,
M = 8,
argvals = argvalsSetting1,
argvalsfd,
uniExpansions = uniExpansions,
B){
set.seed(seed)
## MFPCA
MRSE <- numeric(B)
ErrValues <- matrix(NA, nrow = B, ncol = M)
ErrFuns <- matrix(NA, nrow = B, ncol = M)
## fda
MRSE_RS <- numeric(B)
ErrValuesRS <- matrix(NA, nrow = B, ncol = M)
ErrFunsRS <- matrix(NA, nrow = B, ncol = M)
for (i in 1:B) {
SimData <- simMultiFunData(type = "split",
argvals = argvalsSetting1,
M = M,
eFunType = "Fourier",
eValType = eValType,
N = N)
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
if (s2 > 0) {
SimData2 <- addError(SimData$simData, sd = sqrt(s2))
} else {
SimData2 <- SimData$simData
}
MFPCAObj <- MFPCA(SimData2,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(SimData2 - fits, squared = TRUE)
den <- norm(SimData2, squared = TRUE)
MRSE[i] <- mean(num / den)
# Eigenvalue relative error
ErrValues[i, ] <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns[i, ] <- pmin(ErrFunsA, ErrFunsB)
### fda package
fdData1 <- aperm(SimData2@.Data[[1]]@X, c(2, 1))
fdData2 <- aperm(SimData2@.Data[[2]]@X, c(2, 1))
fdData <- array(NA, dim = c(nrow(fdData1), ncol(fdData1), 2))
fdData[ , , 1] <- fdData1
fdData[ , , 2] <- fdData2
fdarange <- range(argvalsfd)
fdabasis = create.bspline.basis(rangeval = fdarange,
nbasis = 15,
norder = 4)
fdatime <- argvalsfd
fdafd <- smooth.basis(fdatime, fdData, fdabasis)$fd
nharm <- M
fdapcaList <- pca.fd(fdafd, nharm)
fdaValues <- fdapcaList$values[1:M]
harmEval <- eval.fd(argvalsfd, fdapcaList$harmonics)
fdaFuns <- multiFunData(list(funData(argvals = list(argvalsfd),
X = t(harmEval[, , 1])),
funData(argvals = list(argvalsfd),
X = t(harmEval[, , 2]))))
sgn <- ifelse(scalarProduct(TrueFunctions, fdaFuns) < 0, -1, 1)
for (j in 1:2) {
fdaFuns[[j]]@X <- fdaFuns[[j]]@X * sgn
}
### MRSE
meanEval <- eval.fd(argvalsfd, fdapcaList$meanfd)
scores_RS <- apply(fdapcaList$scores[, , , drop = FALSE], c(1, 2), sum)
fits_fd <- array(0, dim = c(length(argvalsfd), N, 2))
for (j in 1:2) {
for (m in 1:M) {
fits_fd[, , j] <- fits_fd[, , j] + outer(harmEval[, m, j], scores_RS[, m])
}
fits_fd[, , j] <- sweep(fits_fd[, , j], 1, meanEval[, , j], "+")
}
fits_RS <- multiFunData(list(funData(argvals = list(argvalsfd),
X = t(fits_fd[, , 1])),
funData(argvals = list(argvalsfd),
X = t(fits_fd[, , 2]))))
numRS <- norm(SimData2 - fits_RS, squared = TRUE)
denRS <- norm(SimData2, squared = TRUE)
MRSE_RS[i] <- mean(numRS / denRS)
# Eigenvalue relative error
ErrValuesRS[i, ] <- (TrueValues - fdaValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsARS <- norm(TrueFunctions - fdaFuns, squared = TRUE)
ErrFunsBRS <- norm(TrueFunctions + fdaFuns, squared = TRUE)
ErrFunsCRS <- norm(TrueFunctions, squared = TRUE)
ErrFunsRS[i, ] <- pmin(ErrFunsARS, ErrFunsBRS)
}
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns,
MRSERS = MRSE_RS,
ErrValuesRS = ErrValuesRS,
ErrFunsRS = ErrFunsRS))
}A continuación, se define el número de observaciones a simular, el error, el dominio sobre el cual se generan las observaciones y el número de componentes principales univariadas a estimar. Además, se muestra un ejemplo de un dataset generado con valores propios que decaen exponencialmente.
N <- 250
s2 <- 0.25 # measurement error
argvalsSetting1 <- list(seq(0, 1, length = 100), # 100 sampling points
seq(0, 1, length = 100))
argvalsfd <- seq(0, 1, length = 100)
uniExpansions <- list(list(type = "uFPCA", npc = 8),
list(type = "uFPCA", npc = 8))
M <- 8SimDataS1 <- simMultiFunData(type = "split",
argvals = argvalsSetting1,
M = M,
eFunType = "Fourier",
eValType = "exponential",
N = N)
SimDataS1s2 <- addError(SimDataS1$simData, sd = sqrt(s2))
A continuación, se realiza la simulación de esta primera configuración \(B=100\) veces con ambos tipos de decaimiento de los valores propios (exponencial y lineal) y errores de medida (\(\sigma^2=0\) y \(\sigma^2=0.25\)).
S1_exp_s0 <- Sim_Setting1(eValType = "exponential",
s2 = 0,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
argvalsfd = argvalsfd,
uniExpansions = uniExpansions,
B = 100)
S1_exp_s <- Sim_Setting1(eValType = "exponential",
s2 = s2,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
argvalsfd = argvalsfd,
uniExpansions = uniExpansions,
B = 100)
S1_lin_s0 <- Sim_Setting1(eValType = "linear",
s2 = 0,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
argvalsfd = argvalsfd,
uniExpansions = uniExpansions,
B = 100)
S1_lin_s <- Sim_Setting1(eValType = "linear",
s2 = s2,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
uniExpansions = uniExpansions,
argvalsfd = argvalsfd,
B = 100)Los resultados se presentan en la Figura 18.
Para la segunda configuración, se evalúa el rendimiento de la nueva propuesta para datos dispersos con las primeras \(M=8\) funciones en la base de Fourier. Se realiza la comparación sin dispersión (Full data), dispersión media (Medium sparsity) y dispersión alta (High sparsity). A diferencia de la configuración anterior, se consideran realizaciones de un proceso funcional trivariado con distintos dominios (\(\mathcal{T}_1=[-1,-0.5]\), \(\mathcal{T}_2=[0,1]\), _3=[1.5,2]) y número de puntos discretizados distintos (\(S_1=S3=50\) y \(S_2=100\)). Además, para cada opción de sparsity se tienen las siguientes aclaraciones:
El MFPCA se calcula dependiendo de la dispersión. Para las simulaciones de datos completos y dispersión media se usan \(M_1=M_2=M_3=8\) componentes principales funcionales univariadas. Con dispersión alta, se usan \(M_1=M_3=3\) y \(M_2=5\) por estabilidad.
A continuación, se muestra un ejemplo de un dataset simulado con un decaimiento exponencial de los valores propios. Al igual que en el caso anterior, se hace uso de la función simMultiFunData().
N <- 250
s2 <- 0.25 # measurement error
T1 <- seq(-1, -0.5, length = 50)
T2 <- seq(0, 1, length = 100)
T3 <- seq(1.5, 2, length = 50)
M <- 8
Ms <- c(8, 8, 8)
MsHigh <- c(3, 5, 3)
argvalsSetting2 <- list(T1, T2, T3)
M1 <- Ms[1]
M2 <- Ms[2]
M3 <- Ms[3]
M1High <- MsHigh[1]
M2High <- MsHigh[2]
M3High <- MsHigh[3]
uniExpansions <- list(list(type = "uFPCA", npc = M1),
list(type = "uFPCA", npc = M2),
list(type = "uFPCA", npc = M3))
uniExpansionsHigh <- list(list(type = "uFPCA", npc = M1High),
list(type = "uFPCA", npc = M2High),
list(type = "uFPCA", npc = M3High))
SimDataS2 <- simMultiFunData(type = "split",
argvals = argvalsSetting2,
M = M,
eFunType = "Fourier",
eValType = "exponential",
N = N)
SimDataS2s2 <- addError(SimDataS2$simData, sd = sqrt(s2))
Para la simulación, primero se definen unas funciones para cada tipo de dispersión a considerar en la simulación.
FullSparsity <- function(Tlens, SimData, SimData2, M, uniExpansions){
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
TrueObs <- SimData2
SimData3 <- SimData2
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(TrueObs - fits, squared = TRUE)
den <- norm(TrueObs, squared = TRUE)
MRSE <- mean(num / den)
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}
MediumSparsity <- function(Tlens, SimData, SimData2, M, uniExpansions){
TrueObs <- SimData2
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
minObs <- round(0.3 * Tlens)
maxObs <- round(0.5 * Tlens)
SimData3 <- sparsify(SimData2, minObs = minObs, maxObs = maxObs)
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(TrueObs - fits, squared = TRUE)
den <- norm(TrueObs, squared = TRUE)
MRSE <- mean(num / den)
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}
HighSparsity <- function(Tlens, SimData, SimData2, M, uniExpansions, MsHigh){
TrueObs <- SimData2
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
minObs <- pmax(round(0.1 * Tlens), MsHigh + 1)
maxObs <- round(0.3 * Tlens)
SimData3 <- sparsify(SimData2, minObs = minObs, maxObs = maxObs)
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(TrueObs - fits, squared = TRUE)
den <- norm(TrueObs, squared = TRUE)
MRSE <- mean(num / den)
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}La función para realizar las simulaciones se presenta a continuación.
Sim_Setting2 <- function(eValType = c("exponential", "linear"),
s2 = 0,
seed = NULL,
N = 250,
M,
Ms,
MsHigh,
argvals = argvalsSetting2,
B){
M1 <- Ms[1]
M2 <- Ms[2]
M3 <- Ms[3]
M1High <- MsHigh[1]
M2High <- MsHigh[2]
M3High <- MsHigh[3]
uniExpansions <- list(list(type = "uFPCA", npc = M1),
list(type = "uFPCA", npc = M2),
list(type = "uFPCA", npc = M3))
uniExpansionsHigh <- list(list(type = "uFPCA", npc = M1High),
list(type = "uFPCA", npc = M2High),
list(type = "uFPCA", npc = M3High))
Tlens <- unlist(lapply(argvals, length))
MRSE_Full <- numeric(B)
ErrValues_Full <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Full <- matrix(NA, nrow = B, ncol = M)
MRSE_Medium <- numeric(B)
ErrValues_Medium <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Medium <- matrix(NA, nrow = B, ncol = M)
MRSE_High <- numeric(B)
ErrValues_High <- matrix(NA, nrow = B, ncol = M)
ErrFuns_High <- matrix(NA, nrow = B, ncol = M)
set.seed(seed)
## MFPCA
for (i in 1:B) {
SimData <- simMultiFunData(type = "split",
argvals = argvals,
M = M,
eFunType = "Fourier",
eValType = eValType,
N = N)
if (s2 > 0) {
SimData2 <- addError(SimData$simData, sd = sqrt(s2))
} else {
SimData2 <- SimData$simData
}
Full <- FullSparsity(Tlens = Tlens,
SimData = SimData,
SimData2 = SimData2,
M = M,
uniExpansions = uniExpansions)
Medium <- MediumSparsity(Tlens = Tlens,
SimData = SimData,
SimData2 = SimData2,
M = M,
uniExpansions = uniExpansions)
High <- HighSparsity(Tlens = Tlens,
SimData = SimData,
SimData2 = SimData2,
M = M,
uniExpansions = uniExpansionsHigh,
MsHigh = MsHigh)
MRSE_Full[i] <- Full$MRSE
ErrValues_Full[i, ] <- Full$ErrValues
ErrFuns_Full[i, ] <- Full$ErrFuns
MRSE_Medium[i] <- Medium$MRSE
ErrValues_Medium[i, ] <- Medium$ErrValues
ErrFuns_Medium[i, ] <- Medium$ErrFuns
MRSE_High[i] <- High$MRSE
ErrValues_High[i, ] <- High$ErrValues
ErrFuns_High[i, ] <- High$ErrFuns
}
MRSE = list(Full = MRSE_Full,
Medium = MRSE_Medium,
High = MRSE_High)
ErrValues = list(Full = ErrValues_Full,
Medium = ErrValues_Medium,
High = ErrValues_High)
ErrFuns <- list(Full = ErrFuns_Full,
Medium = ErrFuns_Medium,
High = ErrFuns_High)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}Y se realizan las respectivas simulaciones.
N <- 250
s2 <- 0.25 # measurement error
T1 <- seq(-1, -0.5, length = 50)
T2 <- seq(0, 1, length = 100)
T3 <- seq(1.5, 2, length = 50)
M <- 8
Ms <- c(8, 8, 8)
MsHigh <- c(3, 5, 3)
argvalsSetting2 <- list(T1, T2, T3)
S2_exp_s0 <- Sim_Setting2(eValType = "exponential",
s2 = 0,
seed = 2,
N = N,
M = M,
Ms = Ms,
MsHigh = MsHigh,
argvals = argvalsSetting2,
B = 100)
S2_exp_s <- Sim_Setting2(eValType = "exponential",
s2 = s2,
seed = 2,
N = 250,
M = 8,
Ms = rep(8, 3),
MsHigh = c(3, 5, 3),
argvals = argvalsSetting2,
B = 100)
S2_lin_s0 <- Sim_Setting2(eValType = "linear",
s2 = 0,
seed = 2,
N = 250,
M = 8,
Ms = rep(8, 3),
MsHigh = c(3, 5, 3),
argvals = argvalsSetting2,
B = 100)
S2_lin_s <- Sim_Setting2(eValType = "linear",
s2 = s2,
seed = 2,
N = 250,
M = 8,
Ms = rep(8, 3),
MsHigh = c(3, 5, 3),
argvals = argvalsSetting2,
B = 100)Los resultados se presentan en la Figura 22.
Expansions <- function(SimData, SimData2, TrueFuns, M, uniExpansions){
TrueValues <- SimData$trueVals[1:M]
TrueFunctions <- TrueFuns
SimData3 <- SimData2
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(ErrValues = ErrValues,
ErrFuns = ErrFuns))
}Sim_Setting3 <- function(eValType = "exponential",
s2 = 0,
seed = NULL,
N = 250,
M,
Ms,
argvals = argvalsSetting3,
B){
uniExpansionsBasis <- list(list(type = "splines2D",
k = c(10, 12)),
list(type = "splines1D",
k = 15))
if(s2 > 0){
uniExpansionsBasis <- list(list(type = "splines2D",
k = c(10, 12)),
list(type = "splines1Dpen",
k = 15))
}
alphaRange = list(v = c(10^-5, 10^5),
w = c(10^-5, 10^5))
uniExpansionsTensor <- list(list(type = "FCP_TPA",
npc = 20,
alphaRange = alphaRange),
list(type = "uFPCA",
npc = 15))
ErrValues_Splines <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Splines <- matrix(NA, nrow = B, ncol = M)
ErrValues_Tensor <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Tensor <- matrix(NA, nrow = B, ncol = M)
set.seed(seed)
## MFPCA
pb <- txtProgressBar(min = 0, max = B, style = 3)
on.exit(close(pb), add = TRUE)
for (i in 1:B) {
SimData <- simMultiFunData(type = "weighted",
argvals = argvals,
M = Ms,
eFunType = list(c("Fourier", "Fourier"),
c("Poly")),
eValType = eValType,
N = N)
TrueFuns <- multiFunData(
funData(argvals = SimData$trueFuns[[1]]@argvals,
X = SimData$trueFuns[[1]]@X[1:M, , , drop = FALSE]),
funData(argvals = SimData$trueFuns[[2]]@argvals,
X = SimData$trueFuns[[2]]@X[1:M, , drop = FALSE])
)
if (s2 > 0) {
SimData2 <- addError(SimData$simData, sd = sqrt(s2))
} else {
SimData2 <- SimData$simData
}
Splines <- Expansions(SimData = SimData,
SimData2 = SimData2,
TrueFuns = TrueFuns,
M = M,
uniExpansions = uniExpansionsBasis)
Tensor <- Expansions(SimData = SimData,
SimData2 = SimData2,
TrueFuns = TrueFuns,
M = M,
uniExpansions = uniExpansionsTensor)
ErrValues_Splines[i, ] <- Splines$ErrValues
ErrFuns_Splines[i, ] <- Splines$ErrFuns
ErrValues_Tensor[i, ] <- Tensor$ErrValues
ErrFuns_Tensor[i, ] <- Tensor$ErrFuns
setTxtProgressBar(pb, i)
}
ErrValues = list(Splines = ErrValues_Splines,
Tensor = ErrValues_Tensor)
ErrFuns <- list(Splines = ErrFuns_Splines,
Tensor = ErrFuns_Tensor)
return(list(ErrValues = ErrValues,
ErrFuns = ErrFuns))
}T1_1 <- seq(0, 1, length = 100)
T1_2 <- seq(0, 0.5, length = 50)
T2 <- seq(-1, 1, length = 200)
argvalsSetting3 <- list(list(T1_1, T1_2),
list(T2))
Ms <- list(c(5, 5), 25)
M <- 12
N <- 250
s2 <- 0.25
S3_exp_s0 <- Sim_Setting3(eValType = "exponential",
s2 = 0,
seed = 3,
N = 250,
M,
Ms,
argvals = argvalsSetting3,
B = 100)
S3_exp_s <- Sim_Setting3(eValType = "exponential",
s2 = s2,
seed = 3,
N = 250,
M,
Ms,
argvals = argvalsSetting3,
B = 100)