CovFIntegral {frechet}

Generalized Fréchet integrals of covariance matrix

Description

Calculating generalized Fréchet integrals of covariance (equipped with Frobenius norm)

Usage

CovFIntegral(phi, t_out, X)

Arguments

phi

An eigenfunction along which we want to project the network

t_out

Support of phi

X

A three dimensional array of dimension length(t_out) x m x m, where X[i,,] is an m x m covariance matrix.

Value

A list of the following:

f

An adjacency matrix which corresponds to the Fréchet integral of X along phi

References

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.

Examples

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)