In statKim/robfpca: Functional principal component analysis for partially observed elliptical process

knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
library(robfpca)

Functional PCA for partially observed elliptical process

This is the R package robfpca implementing the robust functional principal component analysis (FPCA) for partially observed functional data from the following paper:

Yeonjoo Park, Hyunsung Kim, and Yaeji Lim (2023). Functional principal component analysis for partially observed elliptical process, Computational Statistics & Data Analysis, 184, 107745. https://doi.org/10.1016/j.csda.2023.107745

The proposed robust FPCA method implements FPCA based on the conditional expectation for the robust covariance function estimate. The robust covariance function is obtained via robust pairwise computation based on Orthogonalized Gnanadesikan-Kettenring (OGK) estimation.

Installation

# install.packages("devtools")
devtools::install_github("statKim/robfpca")

Example

Generate partially observed functional data from heavy-tailed distribution

First, we generate 100 partially observed curves from heavy-tailed $t_3$ distribution.

library(robfpca)

# Generate partially observed curves from heavy-tailed distribution
set.seed(46)
X.list <- sim_delaigle(n = 100,
                       type = "partial",
                       dist = "tdist")
X <- list2matrix(X.list)   # transform list to matrix
gr <- seq(0, 1, length.out = ncol(X))   # observed timepoints
matplot(gr, t(X), 
        type = "l",
        xlab = "t", ylab = "X(t)")

Robust FPCA for partially observed functional data

# Proposed FPCA with bandwidth = 0.3
robfpca.obj <- robfpca.partial(X,
                               type = "huber",
                               # PVE = 0.99,
                               K = 4,   # we use true dimension
                               bw = 0.3)
fpc.score <- robfpca.obj$pc.score  # FPC scores

# First three eigenfunctions
eig.true <- get_delaigle_eigen(gr, model = 2)   # True eigenfunctions
eig.robfpca <- check_eigen_sign(robfpca.obj$eig.fun, 
                                eig.true)   # match eigen directions
par(mfrow = c(1, 3))
for (i in 1:3) {
    plot(gr, eig.true[, i],
         type = "l",
         lwd = 2,
         ylim = c(-2, 2.5),
         xlab = "t", 
         ylab = paste("PC", i),
         main = paste("Eigenfunction", i))
    lines(gr, eig.robfpca[, i], col = 2, lwd = 2)

    if (i == 1) {
        legend("bottomright", 
               c("True","Proposed"),
               lty = c(1, 1),
               col = c(1, 2))
    }
}

Predict FPC score

# Select 4 curves that have sufficiently long missing periods
set.seed(46)
idx <- sample(which(rowSums(is.na(X)) > 10), 4)
new_data <- X[idx, ]   # example of new data

# new_data <- X[c(1,4,5), ]   # example of new data

# Predict the FPC scores
pred_score <- predict(robfpca.obj, type = "score", newdata = new_data)
pred_score

Reconstruction

pred_reconstr <- predict(robfpca.obj, type = "reconstr", newdata = new_data)
pred_reconstr[1, ]  # reconstructed curve of 1st new_data

par(mfrow = c(1, 2))
matplot(t(new_data), type = "l",
        ylim = range(pred_reconstr, new_data, na.rm = T) + c(-0.5, 0.5),
        xlab = "t", ylab = "", main = "Observed curves")
matplot(t(pred_reconstr), type = "l",
        ylim = range(pred_reconstr, new_data, na.rm = T) + c(-0.5, 0.5),
        xlab = "t", ylab = "", main = "Reconstructed curves")

Completion

pred_comp <- predict(robfpca.obj, type = "comp", newdata = new_data)
pred_comp[1, ]  # predicted missing parts of 1st new_data

matplot(t(new_data), type = "l",
        ylim = range(pred_comp, new_data, na.rm = T) + c(-0.5, 0.5),
        xlab = "t", ylab = "", main = "Completion")
matlines(t(pred_comp), type = "l", lty = 1, lwd = 2)

Reproducible Simulation Results

If you use simulation.R in Here, you can obtain reproducible results in our paper.

Additional example

In example, you can also check the comparison result between Proposed FPCA and Sparse FPCA.

statKim/robfpca documentation built on April 15, 2023, 10:12 p.m.