# R/covZ.r In face: Fast Covariance Estimation for Sparse Functional Data

covZ <- function(t_seq, Y, C, sigma2)
{
#########################################################################
#     covZ: calculate the covariance matrix of paired-y product
#
#     Input
#     t_seq: defined grid point m-by-1
#     Y: response vector m-by-1
#     C: covariance matrix of X m-by-m
#     sigma2: variance of random error
#
#     Output
#     CovZZ: covariance matrix of paired-y product m^2-by-m^2
#     Z.info: vectorized Z.tri and  the corresponding grid points
#     Z.tri: upper triangle matrix of paired-y product m-by-m
#
#########################################################################
temp <- expand.grid(Y,Y)
t.pair <- expand.grid(t_seq,t_seq)
d <- length(Y)
id <- expand.grid(1:d,1:d)
m <- d

Z <- apply(temp, 1, function(z) prod(z))
Z.info <- as.matrix(cbind(Z,t.pair,id))
Z.info <- Z.info[Z.info[,5] <= Z.info[,4],]

#   Z.mat <- matrix(Z,nrow=length(Y),ncol=length(Y),byrow=T)
#   Z.tri <- Z.mat
#   Z.tri[lower.tri(Z.tri)] <- NA

Z.tri <- matrix(NA,nrow=length(Y),ncol=length(Y))
Z.tri[lower.tri(Z.tri, diag=T)] <- Z.info[,1]
Z.tri <- t(Z.tri)

# Quick-n-dirty implementation (for testing only)
#   CovZZ <- matrix(NA,nrow=length(Z),ncol=length(Z))
#   time.start <- proc.time()
#   for (i in 1:m)
#   {
#     for (j in 1:m)
#     {
#       row.id <- d * (i - 1) + j
#       for (k in 1:m)
#       {
#         for (l in 1:m)
#         {
#           col.id <- d * (k - 1) + l
#           if (row.id <= col.id)
#           {
#             CovZZ[row.id,col.id] = C[i,k] * C[j,l] + C[i,l] * C[j,k] +
#               (i==k) * (j==l) *sigma2^2 + (i==l) * (j==k) *sigma2^2 +
#               (C[i,k] * (j==l) + C[i,l] * (j==k) + C[j,k] * (i==l)
#                + C[j,l] * (i==k)) * sigma2
#           }
#
#         }
#       }
#     }
#   }
#   print((proc.time()-time.start)[3])
#
#   # Avoid loops (even slower, for testing only)
#   time.start <- proc.time()
#   id1 <- as.matrix(expand.grid(1:d,1:d,1:d,1:d))
#   id2 <- as.matrix(expand.grid(1:d^2,1:d^2))
#   id.tab <- as.matrix(cbind(id1,id2))
#   id.tab <- id.tab[id.tab[,5] <= id.tab[,6],]
#   CovZZ2.vec = apply(id.tab,1,function(x)
#     C[x[2],x[4]] * C[x[1],x[3]] + C[x[2],x[3]] * C[x[1],x[4]] +
#     (x[2]==x[4]) * (x[1]==x[3]) *sigma2^2 +
#     (x[2]==x[3]) * (x[1]==x[4]) *sigma2^2 +
#     (C[x[2],x[4]] * (x[1]==x[3]) + C[x[2],x[3]] * (x[1]==x[4]) +
#        C[x[1],x[4]] * (x[2]==x[3]) + C[x[1],x[3]] * (x[2]==x[4])) * sigma2)
#
#   CovZZ2 <- matrix(NA,nrow=length(Z),ncol=length(Z))
#   CovZZ2[upper.tri(CovZZ2, diag=T)] <- CovZZ2.vec
#   print((proc.time()-time.start)[3])

# Alternative vectorizing
#time.start <- proc.time()
id1 <- as.matrix(expand.grid(1:d,1:d,1:d,1:d))
id2 <- as.matrix(expand.grid(1:d^2,1:d^2))
id.tab <- as.matrix(cbind(id1,id2))
id.tab <- id.tab[id.tab[,5] <= id.tab[,6],] # symmetric
# kronecker delta terms
# \delta_{jk},\delta_{j'k'},\delta_{jk'},\delta_{j'k}
k.del <- cbind(id.tab[,1]==id.tab[,3],id.tab[,2]==id.tab[,4],
id.tab[,1]==id.tab[,4],id.tab[,2]==id.tab[,3])
# Cov operator terms
# C_{jk},C_{j'k'},C_{jk'},C_{j'k}
c.term <- cbind(C[cbind(id.tab[,1],id.tab[,3])],
C[cbind(id.tab[,2],id.tab[,4])],
C[cbind(id.tab[,1],id.tab[,4])],
C[cbind(id.tab[,2],id.tab[,3])])
# Form a matrix with pre-calculated terms
id.tab <- cbind(id.tab,k.del,c.term)
# CovZZ = C_{jk}C_{j'k'} + C_{jk'}C_{j'k}
# + (\delta_{jk}\delta_{j'k'} + \delta_{jk'}\delta_{j'k})\sigma_2^2 +
# (C_{jk}\delta_{j'k'} + C_{jk'}\delta_{j'k} + C_{j'k}\delta_{jk'} +
# C_{j'k'}\delta_{jk})\sigma_2
CovZZ3.vec = id.tab[,12] * id.tab[,11] + id.tab[,14] * id.tab[,13] +
(id.tab[,8] * id.tab[,7] + id.tab[,10] * id.tab[,9]) *sigma2^2 +
(id.tab[,12] * id.tab[,7] + id.tab[,14] * id.tab[,9] +
id.tab[,13] * id.tab[,10] + id.tab[,11] * id.tab[,8]) * sigma2
# CovZZ.vec = vech(CovZZ)
CovZZ3 <- matrix(NA,nrow=length(Z),ncol=length(Z))
CovZZ3[upper.tri(CovZZ3, diag=T)] <- CovZZ3.vec
#print((proc.time()-time.start)[3])

return (list(CovZZ = CovZZ3,CovZZ.vech = CovZZ3.vec,
Z.info = Z.info,Z.tri = Z.tri))

}


## Try the face package in your browser

Any scripts or data that you put into this service are public.

face documentation built on May 2, 2019, 6:47 a.m.