R/CovFIntegral.R

In frechet: Statistical Analysis for Random Objects and Non-Euclidean Data

#' @title Generalized Fréchet integrals of covariance matrix
#' @description Calculating generalized Fréchet integrals of covariance (equipped with Frobenius norm) 
#' @param phi An eigenfunction along which we want to project the network
#' @param t_out Support of \code{phi}
#' @param X A three dimensional array of dimension \code{length(t_out) x m x m}, where \code{X[i,,]} is an m x m covariance matrix. 
#' @return A list of the following:
#' \item{f}{An adjacency matrix which corresponds to the Fréchet integral of \code{X} along \code{phi}}
#' @examples
#' \donttest{
#' set.seed(5)
#' library(mpoly)
#' n <- 100
#' N <- 50
#' t_out <- seq(0,1,length.out = N)
#' 
#' p2 <- as.function(mpoly::jacobi(2,4,3),silent=TRUE)
#' p4 <- as.function(mpoly::jacobi(4,4,3),silent=TRUE)
#' p6 <- as.function(mpoly::jacobi(6,4,3),silent=TRUE)
#' 
#' # first three eigenfunctions
#' phi1 <- function(t){
#'   p2(2*t-1)*t^(1.5)*(1-t)^2/(integrate(function(x) p2(2*x-1)^2*x^(3)*(1-x)^4,0,1))$value^(1/2)
#' }
#' phi2 <- function(t){
#'   p4(2*t-1)*t^(1.5)*(1-t)^2/(integrate(function(x) p4(2*x-1)^2*x^(3)*(1-x)^4,0,1))$value^(1/2)
#' }
#' phi3 <- function(t){
#'   p6(2*t-1)*t^(1.5)*(1-t)^2/(integrate(function(x) p6(2*x-1)^2*x^(3)*(1-x)^4,0,1))$value^(1/2)
#' }
#' 
#' # random component of covariance matrices
#' P12 <- 0 ## elements between communities
#' Score <- matrix(runif(n*4), nrow = n)
#' # first community
#' P1_cov <- 0.5 + 0.4*Score[,1] %*% t(phi1(t_out)) + 0.1*Score[,2] %*% t(phi3(t_out))  
#' # second community
#' P2_cov <- 0.5 + 0.3*Score[,3] %*% t(phi2(t_out)) + 0.1*Score[,4] %*% t(phi3(t_out))  
#' P1_diag <- 2 #diagonal elements of the first community
#' P2_diag <- 3 #diagonal elements of the second community
#' 
#' # create Network edge matrix 
#' N_net1 <- 5 # first community number
#' N_net2 <- 5 # second community number
#' 
#' # I: four dimension array of n x n matrix of squared distances between the time point u  
#' # of the ith process and process and the time point v of the jth object process, 
#' # e.g.: I[i,j,u,v] <- d_F^2(X_i(u) X_j(v)).
#' I <- array(0, dim = c(n,n,N,N))
#' for(u in 1:N){
#'   for(v in 1:N){
#'     #frobenius norm between two adjcent matrix  
#'     I[,,u,v] <- outer(P1_cov[,u], P1_cov[,v], function(a1, a2) (a1-a2)^2*(N_net1^2-N_net1)) + 
#'       outer(P2_cov[,u], P2_cov[,v], function(a1, a2) (a1-a2)^2*(N_net2^2-N_net2)) 
#'   }
#' }
#' 
#' # check ObjCov work 
#' Cov_result <- ObjCov(t_out, I, 3, smooth=FALSE)
#' Cov_result$lambda  # 0.266 0.15 0.04
#' 
#' # e.g. subj 2 
#' subj <- 2
#' # X_mat is the network for varying times with X[i,,] 
#' # is the adjacency matrices for the ith time point
#' X_mat <- array(0, c(N,(N_net1+N_net2), (N_net1+N_net2)))
#' for(i in 1:N){
#'   # edge between communities is P12
#'   Mat <- matrix(P12, nrow = (N_net1+N_net2), ncol = (N_net1+N_net2)) 
#'   # edge within the first communitiy is P1
#'   Mat[1:N_net1, 1:N_net1] <- P1_cov[subj, i] 
#'   # edge within the second community is P2
#'   Mat[(N_net1+1):(N_net1+N_net2), (N_net1+1):(N_net1+N_net2)] <- P2_cov[subj, i] 
#'   diag(Mat) <- c(rep(P1_diag,N_net1),rep(P2_diag, N_net2)) #diagonal element is 0 
#'   X_mat[i,,] <- Mat
#' }
#' # output the functional principal network(adjacency matrice) of the second eigenfunction
#' CovFIntegral(Cov_result$phi[,2], t_out, X_mat)
#' }
#' @references 
#' \cite{Dubey, P., & Müller, H. G. (2020). Functional models for time‐varying random objects. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(2), 275-327.}
#' @import fdapace
#' @importFrom Matrix nearPD
#' @export 

CovFIntegral <- function(phi, t_out, X){
  
  N <- length(phi)
  m <- dim(X)[2]
  if(dim(X)[1] != N){
    stop("length of first argument of X and length of phi are not consistenet ")
  }
  if(dim(X)[3] != dim(X)[2]){
    stop("X[i,,] should be a square matrix")
  }
  
  phi_out <- phi / fdapace::trapzRcpp(t_out, phi)
  g_mini <- matrix(0, nrow = m, ncol = m)
  for(i in 1:(m-1)){
    for(j in (i+1):m){
      g_mini[i,j] <- fdapace::trapzRcpp(t_out, phi_out * X[,i,j])
      g_mini[j,i] <- g_mini[i,j]
    }
    g_mini[i,i] <- max(fdapace::trapzRcpp(t_out, phi * X[,i,i]),10^(-5))
  }
  g_mini[m,m] <- max(fdapace::trapzRcpp(t_out, phi * X[,m,m]),10^(-5))
  # find the nearest positive definite matrix to g_mini
  f = Matrix::nearPD(g_mini,doSym = TRUE)$mat
  return(list(f=f))
}