Nothing
## -----------------------------------------------------------------------------
library(Tlasso)
## -----------------------------------------------------------------------------
m.vec = c(5,5,5) # dimensionality of a tensor
# m1, m2, m3
n = 5 # sample size
Omega.true.list = list()
for (k in 1:length(m.vec)) {
Omega.true.list[[k]] = ChainOmega(m.vec[k], sd = k)
}
Omega.true.list[[1]]
## -----------------------------------------------------------------------------
Sigma.true.list = list()
for (k in 1:length(m.vec)) {
Sigma.true.list[[k]] = solve(Omega.true.list[[k]])
} # generate covariance matrices list
DATA=Trnorm(n,m.vec,Sigma.list=Sigma.true.list)
# obersavations from tensor normal distribution
## ---- eval=FALSE--------------------------------------------------------------
# DATA2=Trnorm(n,m.vec)
# # default is triangle graph
# # equivalent to DATA2 = Trnorm(n,m.vec, type='Chain', sd=1)
# DATA3=Trnorm(n,m.vec,type='Neighbor')
# # 4 nearest-neighbor graph
# # equivalent to DATA3 = Trnorm(n,m.vec, type='Neighbor', sd=1, knn=4)
## -----------------------------------------------------------------------------
lambda.thm = 20*c( sqrt(log(m.vec[1])/(n*prod(m.vec))),
sqrt(log(m.vec[2])/(n*prod(m.vec))),
sqrt(log(m.vec[3])/(n*prod(m.vec))))
# lambda.thm is regularization parameter
out.tlasso = Tlasso.fit(DATA,T=1,lambda.vec = lambda.thm)
# output is a list of estimation of precision matrices
out.tlasso[[1]]
## -----------------------------------------------------------------------------
# compare out.tlasso and Omega.true.list
# main diagnoal is taken into consideration
est.analysis(out.tlasso,Omega.true.list,offdiag=FALSE)
## -----------------------------------------------------------------------------
mat.list=list() # list of matrices of test statistic value
for ( k in 1:length(m.vec)) {
rho=covres(DATA, out.tlasso, k = k)
# sample covariance matrix of residuals, including diagnoal
bias_rho=biascor(rho,out.tlasso,k=k)
# bias corrected sample covariance of residuals, excluding diagnoal
varpi2=varcor(DATA, out.tlasso, k = k)
# variance correction term for kth mode's sample covariance of residuals
tautest=matrix(0,m.vec[k],m.vec[k])
for( i in 1:(m.vec[k]-1)) {
for ( j in (i+1):m.vec[k]){
tautest[j,i]=tautest[i,j]=sqrt((n-1)*prod(m.vec[-k]))*
bias_rho[i,j]/sqrt(varpi2*rho[i,i]*rho[j,j])
# compute final test statistic
}
}
mat.list[[k]]=tautest
}
mat.list[[1]]
## -----------------------------------------------------------------------------
# inference measures (off-diagnoal), critical value is 0.975 quantile of standard normal
infer.analysis(mat.list, qnorm(0.975), Omega.true.list, offdiag=TRUE)
## ---- eval=FALSE--------------------------------------------------------------
# k=1 # interested mode
# upsilon=0.1 # control level
#
# # compute the difference between FDP and upsilon
# fun=function(varsigma,mk,upsilon,tautest) {
# return((2*(1-pnorm(varsigma))*mk*(mk-1))/max(1,sum(sign(abs(tautest) > varsigma))) - upsilon)
# }
#
# # select a critical value in (0,6) that has the samllest difference
# diff=c();ind=1;inter=seq(0,6,0.0001)
# for (varsigma in inter) {
# diff[ind]=fun(varsigma,mk=m.vec[k],upsilon=upsilon,tautest=mat.list[[k]])
# ind=ind+1
# }
# # the smallest critical value that constrains FDP under upsilon
# critical=inter[min(which(diff < 0))]
#
# # testing hypothesis with the critcal value
# # FDR will converge to the limit proved in Lyu et al. 2019.
# inference.FDR=infer.analysis(mat.list, critical, Omega.true.list, offdiag=TRUE)
## ---- fig.height=5 , fig.width=5, fig.align='center'--------------------------
k=1 # interested mode
# true graph structure.
# set thres=0 in case true edge is eliminated
graph.pattern(Omega.true.list[[1]],main='True graph of mode 1',thres=0)
inf.mat=mat.list[[k]] > qnorm(0.975)
# set thres=0 (<1) since inf.mat is logical
graph.pattern(inf.mat,main='Inference of mode 1',thres=0)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.