rm(list=ls())
setwd("~/Dropbox/JournalSubmission/0Stats20230228/rcode/empanalysis/")
library(corpcor)
library(optimx)
library(numDeriv)
newick2phylog<-function (x.tre, call = match.call()) {
complete <- function(x.tre) {
if (length(x.tre) > 1) {
w <- ""
for (i in 1:length(x.tre)) w <- paste(w, x.tre[i],
sep = "")
x.tre <- w
}
ndroite <- nchar(gsub("[^)]", "", x.tre))
ngauche <- nchar(gsub("[^(]", "", x.tre))
if (ndroite != ngauche)
stop(paste(ngauche, "( versus", ndroite, ")"))
if (regexpr(";", x.tre) == -1)
stop("';' not found")
i <- 0
kint <- 0
kext <- 0
arret <- FALSE
if (regexpr("\\[", x.tre) != -1) {
x.tre <- gsub("\\[[^\\[]*\\]", "", x.tre)
}
x.tre <- gsub(" ", "", x.tre)
while (!arret) {
i <- i + 1
if (substr(x.tre, i, i) == ";")
arret <- TRUE
if (substr(x.tre, i, i + 1) == "(,") {
kext <- kext + 1
add <- paste("Ext", kext, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == ",,") {
kext <- kext + 1
add <- paste("Ext", kext, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == ",)") {
kext <- kext + 1
add <- paste("Ext", kext, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == "(:") {
kext <- kext + 1
add <- paste("Ext", kext, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == ",:") {
kext <- kext + 1
add <- paste("Ext", kext, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == "),") {
kint <- kint + 1
add <- paste("I", kint, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == "))") {
kint <- kint + 1
add <- paste("I", kint, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == "):") {
kint <- kint + 1
add <- paste("I", kint, sep = "")
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
else if (substr(x.tre, i, i + 1) == ");") {
add <- "Root"
x.tre <- paste(substring(x.tre, 1, i), add, substring(x.tre,
i + 1), sep = "")
i <- i + 1
}
}
lab.points <- strsplit(x.tre, "[(),;]")[[1]]
lab.points <- lab.points[lab.points != ""]
no.long <- (regexpr(":", lab.points) == -1)
if (all(no.long)) {
lab.points <- paste(lab.points, ":", c(rep("1", length(no.long) -
1), "0.0"), sep = "")
}
else if (no.long[length(no.long)]) {
lab.points[length(lab.points)] <- paste(lab.points[length(lab.points)],
":0.0", sep = "")
}
else if (any(no.long)) {
stop("Non convenient ancdes.ancdes.array leaves or nodes with and without length")
}
w <- strsplit(x.tre, "[(),;]")[[1]]
w <- w[w != ""]
leurre <- make.names(w, unique = TRUE)
leurre <- gsub("[.]", "_", leurre)
for (i in 1:length(w)) {
old <- paste(w[i])
x.tre <- sub(old, leurre[i], x.tre, fixed = TRUE)
}
w <- strsplit(lab.points, ":")
label <- function(x) {
lab <- x[1]
lab <- gsub("[.]", "_", lab)
return(lab)
}
longueur <- function(x) {
long <- x[2]
return(long)
}
labels <- unlist(lapply(w, label))
longueurs <- unlist(lapply(w, longueur))
labels <- make.names(labels, TRUE)
labels <- gsub("[.]", "_", labels)
w <- labels
for (i in 1:length(w)) {
new <- w[i]
x.tre <- sub(leurre[i], new, x.tre)
}
cat <- rep("", length(w))
for (i in 1:length(w)) {
new <- w[i]
if (regexpr(paste("\\)", new, sep = ""), x.tre) !=
-1)
cat[i] <- "int"
else if (regexpr(paste(",", new, sep = ""), x.tre) !=
-1)
cat[i] <- "ext"
else if (regexpr(paste("\\(", new, sep = ""), x.tre) !=
-1)
cat[i] <- "ext"
else cat[i] <- "unknown"
}
return(list(tre = x.tre, noms = labels, poi = as.numeric(longueurs),
cat = cat))
}
res <- complete(x.tre)
poi <- res$poi
nam <- res$noms
names(poi) <- nam
cat <- res$cat
res <- list(tre = res$tre)
res$leaves <- poi[cat == "ext"]
names(res$leaves) <- nam[cat == "ext"]
res$nodes <- poi[cat == "int"]
names(res$nodes) <- nam[cat == "int"]
listclass <- list()
dnext <- c(names(res$leaves), names(res$nodes))
listpath <- as.list(dnext)
names(listpath) <- dnext
x.tre <- res$tre
while (regexpr("[(]", x.tre) != -1) {
a <- regexpr("\\([^\\(\\)]*\\)", x.tre)
n1 <- a[1] + 1
n2 <- n1 - 3 + attr(a, "match.length")
chasans <- substring(x.tre, n1, n2)
chaavec <- paste("\\(", chasans, "\\)", sep = "")
nam <- unlist(strsplit(chasans, ","))
w1 <- strsplit(x.tre, chaavec)[[1]][2]
parent <- unlist(strsplit(w1, "[,\\);]"))[1]
listclass[[parent]] <- nam
x.tre <- gsub(chaavec, "", x.tre)
w2 <- which(unlist(lapply(listpath, function(x) any(x[1] ==
nam))))
for (i in w2) {
listpath[[i]] <- c(parent, listpath[[i]])
}
}
res$parts <- listclass
res$paths <- listpath
dnext <- c(res$leaves, res$nodes)
names(dnext) <- c(names(res$leaves), names(res$nodes))
res$droot <- unlist(lapply(res$paths, function(x) sum(dnext[x])))
res$call <- call
class(res) <- "phylog"
#if (!add.tools)
return(res)
#return(newick2phylog.addtools(res))
}
comb<-function(n,r){
return(factorial(n)/(factorial(n-r)*factorial(r)))
}
pair_fcn<-function(tmp){ # return pair for "tmp" sequences.
numl=comb(length(tmp),2)
count=0
posit<-array(0,c(numl))
for(i in 1:length(tmp)){
for(j in 1:length(tmp)){
if(i<j){
count=count+1
posit[count]=paste(c(tmp[i]),c(tmp[j]),sep=",")
}
}
}
return(posit)
}# end of pair_fcn.
pair_array<-function(tmp){#generate pair in array format.
pair <-pair_fcn(tmp)
p_arr<-matrix(c(0),nrow=length(pair),ncol=2)
for(i in 1:length(pair)){
p_arr[i,1]=unlist(strsplit(pair[i],","))[1]
p_arr[i,2]=unlist(strsplit(pair[i],","))[2]
}
return(p_arr)
}# end of function pair_array.
ord_fcn<-function(res,tipnames){ #this function gets the acenstor-descendants relationship.
bmtp<-matrix(rev(res$nde),ncol=1)
rvpt <-rev((res$parts))
rept<-array(0,c(length(rvpt),2))
for(i in 1:length(rvpt)){
rept[i,]=unlist(rvpt[i])
}
cmb<-cbind(bmtp,rept)
brnlen<-res$droot[(length(tipnames)+1):length(res$droot)]
root<-matrix(cmb[1,],nrow=1)
cmb<-cmb[-1,]
brnlen<-brnlen[1:(length(brnlen)-1)]
new_ord<-order(brnlen,decreasing=TRUE)
cmb<-cmb[new_ord,]
cmb<-rbind(root,cmb)
return(cmb)
}# end of function ord_fcn.
getntn<-function(res){# this function gets rid of unnecessarily "_" symbol.
size<-length(res$parts)
relarr<-array(0,c(size,3))
rvpt <-(res$parts)
rept<-array(0,c(length(rvpt),2))
for(i in 1:length(rvpt)){
rept[i,]=unlist(rvpt[i])
}
for(i in 1:size){
relarr[i,1]<-names(res$parts)[i]
}
relarr[,2:3]<-rept
temp<-matrix(0,row<-size)
for(j in 2:3){
for (i in 1: size){
stmp<-unlist(strsplit(relarr[,j][i], "_" ))
temp[i]<-stmp[1]
}
relarr[,j]<-temp
}
ndlen<- res$droot[!(res$droot==max(res$droot))]
nam<-names(ndlen)
ck1<-array(0,c(length(nam)))
count<-0
for (ele in c(nam)){
count<-count+1
len <- length( unlist(strsplit(ele ,"_" )))
if( len==2 ){ck1[count]<-1}
}
ndlen<-ndlen[!ck1]
new_ord<-order(ndlen)
relarr<-relarr[new_ord,]
return(relarr)
}# end of function getntn.
getbrnlen<-function(res){#this function is used to obtain branch length.
ndlen<- res$droot[!(res$droot==1)]
nam<-names(ndlen)
ck1<-array(0,c(length(nam)))
count<-0
for (ele in c(nam)){
count<-count+1
len <- length( unlist(strsplit(ele ,"_" )))
if( len==2 ){ck1[count]<-1}
}
ndlen<-ndlen[!ck1]
ndlen<-sort(ndlen)
ck2<-array(0,c(length(ndlen)))
for(i in 1:(length(ndlen)-1)){
if(abs(ndlen[i]-ndlen[i+1])<10^(-5)){ck2[i]=1}
}
ndlen<-ndlen[!ck2]
brnlen<-array(0,c(length(ndlen)))
tmplen<-ndlen
for(i in 1:(length(brnlen)-1)){
brnlen[i]<-tmplen[i+1]-tmplen[i]
}
brnlen[length(brnlen)] <- 1 -tmplen[(length(tmplen))]
return(brnlen)
}# end of function getbrnlen.
cov_mtx<-function(x,branchlength=branchlength,ancdes.array=ancdes.array,nleaves=nleaves,tipnames=tipnames,model.Index=model.Index){ #now it is a function of bt, h, sigma^2 and sigma_H^2
bt<-x[1]
h<-x[2]
sigma_sq<-x[3]
sigma.H_sq<-x[4]
# bt<-1
h<-0.5
# sigma_sq<-0.68
# sigma.H_sq<-0.68
if(model.Index==1){bt<-1}
if(model.Index==2){sigma.H_sq<-0}
if(model.Index==3){bt<-1;sigma.H_sq<-0}
ins_fcn<-function(ist,sqc){#finds position to insert between two parents, for hybrdization only.
ist<-as.numeric(unlist(strsplit(ist,"X"))[2])
arr<-array(0,c(length(otmp)))
for(i in 1:length(arr)){
arr[i]<-as.numeric(unlist(strsplit(sqc[i],"X"))[2])
}
insp<-which(arr==(ist-1))+1
return(insp)
}
var_fcn<-function(){#return the variance.
for(i in 1:length(otmp)){#use to fill other diagonal.
newi<-which(rownames(mtx)%in%otmp[i])
oldi<-which(rownames(omtx)%in%otmp[i])
mtx[newi,newi]<-omtx[oldi,oldi]
}#fill in old value from omtx exclude the new hyd.
prn1<-tmp[which(tmp%in%ins)-1]#grab elements.
prn2<-tmp[which(tmp%in%ins)+1]
prn1<-which(rownames(omtx) %in% prn1)#grab position according to prn1.
prn2<-which(rownames(omtx) %in% prn2)
vhii<- bt^2*h^2*omtx[prn1,prn1]+bt^2*(1-h)^2*omtx[prn2,prn2]+2*bt^2*h*(1-h)*omtx[prn1,prn2]
hii<-which(!(tmp %in% otmp))#use to insert variance for hyd.
mtx[hii,hii]<-vhii #fill in the diagonal hyd.
return(mtx)
}#formula for insertion hyd variance.
fillspcmtx<-function(){#fill matrix due to sepciation.
elm<-function(){ #use to cut one row of the pair array which the speciation happens.
ck<-c(tmp[nsi],tmp[nsj])
for(i in 1:dim(pn_arr)[1]){
if(sum(pn_arr[i,]==ck)==2){break}
}
return(i)}
pn_arr<-pair_array(tmp)
po_arr<-pair_array(otmp)
#search new speciate position.
nsi<-which(!(tmp %in% otmp))[1]
nsj<-which(!(tmp %in% otmp))[2]
osii<-which(!(otmp %in% tmp))
mtx[nsi,nsj]<- omtx[osii,osii]
#Fill in value: the covariance for 2 speciated species equal the variance of the parent.
pn_arr<-pn_arr[-elm(),]#delete the ancdes.array that is already used.
#The following fills covaraince components by the previous matrix.
while(length(pn_arr[,1])>0){
newi<-which(rownames(mtx) %in% pn_arr[1,1])
newj<-which(rownames(mtx) %in% pn_arr[1,2])
if( tmp[nsi] %in% pn_arr[1,]){
otg<-which(!(pn_arr[1,] %in% tmp[nsi]))
oldi<- which( rownames(omtx) %in% otmp[osii])
oldj<-which(rownames(omtx) %in% pn_arr[1,otg])
}
if( tmp[nsj] %in% pn_arr[1,] ){
otg<-which(!(pn_arr[1,] %in% tmp[nsj]))
oldi<- which( rownames(omtx) %in% otmp[osii])
oldj<-which(rownames(omtx) %in% pn_arr[1,otg])
}
if(!(tmp[nsi] %in% pn_arr[1,]) && !(tmp[nsj] %in% pn_arr[1,])){
#detect common between omtx and mtx.
oldi<-which(rownames(omtx) %in% pn_arr[1,1])
oldj<-which(rownames(omtx) %in% pn_arr[1,2])
}
mtx[newi,newj]<-omtx[oldi,oldj]
pn_arr<-pn_arr[-1,]#delete row.
if(length(pn_arr)==2){pn_arr<-matrix(pn_arr,nrow=1)}
}#end of while loop.
mtx<-mtx+t(mtx)
mtx[nsi,nsi]<-omtx[osii,osii]+ branchlength[length(tmp)-1]
mtx[nsj,nsj]<-omtx[osii,osii]+ branchlength[length(tmp)-1]
dianew<-which(tmp %in% otmp )
diaold<-which(otmp %in% tmp )
for(i in 1:length(dianew)){
mtx[dianew[i],dianew[i]]<-omtx[diaold[i],diaold[i]]+branchlength[length(tmp)-1]
}
return(mtx)
}#end of fillspcmtx.
fillhydmtx<-function(){#fill in value into matrix due to hybridzation.
pn_arr<-pair_array(tmp)
while(length(pn_arr[,1])>0){
newi<-which(rownames(mtx) %in% pn_arr[1,1])
newj<-which(rownames(mtx) %in% pn_arr[1,2])
if (ins %in% pn_arr[1,]){#ins is the hybridized node.
otg<-pn_arr[1,which(!(pn_arr[1,] %in% ins ))]
otgj<-which(rownames(omtx) %in% otg)
#the other guy, could be the hybrdized nodes parent or others.
#find the parent of ins.
prn1<-tmp[which(tmp%in%ins)-1]#grab element.
prn2<-tmp[which(tmp%in%ins)+1]
prn1<-which(rownames(omtx) %in% prn1)#grab position.
prn2<-which(rownames(omtx) %in% prn2)
mtx[newi,newj]<-bt*h*omtx[prn1,otgj] +bt*(1-h)*omtx[prn2,otgj] # cov(X, bt*hX+bt*(1-h)Y) we are going to use h=1/2.
}else{#this is not hyd node, just read from previous mtx.
#only need to read from rownames().
oldi<-which(rownames(omtx) %in% pn_arr[1,1])
oldj<-which(rownames(omtx) %in% pn_arr[1,2])
mtx[newi,newj]<-omtx[oldi,oldj]
}#end of else loop .
pn_arr<-pn_arr[-1,] # delete ancdes.array array after using it.
if(length(pn_arr)==2){pn_arr<-matrix(pn_arr,nrow=1)}
}#end of while loop.
return(mtx)
}#end of fillhydmtx.
#THE MAIN PROGRAM for covariance matrix.
ckins<-FALSE # use to check the hybrdized event.
rept<-ancdes.array[,2:3]# the descedant nodes.
bmtp<-matrix((ancdes.array)[,1],ncol=1) #the acenstor node.
loop<-2
tmp=array(0,c(loop))
if(loop==2){tmp=rept[1,]
otmp<-tmp
mtx<-diag(branchlength[1],c(length(tmp)))
rownames(mtx)<-c(tmp)
colnames(mtx)<-c(tmp)
omtx<-mtx
}#end of loop==2
while(loop<length(bmtp)){#loaded the acenstor-descendant ancdes.array.
loop<-loop+1#use loop to use the ancdes.array
tmp=array(0,c(length(otmp)+1))#the new seq.
mtx<-matrix(0,nrow=length(tmp),ncol=length(tmp))
q=loop-1#index for matching the right element: will use below.
op<-which(otmp==bmtp[q])#index for insertion position.
if(length(op)!=0){#op!=0 means that weve detected speciation.
tmp[op:(op+1)]=rept[q,] #insertion the new speciation species.
if(op==1){tmp[(op+2):length(tmp)]=otmp[(op+1):(length(tmp)-1)]}
if((op+1)==length(tmp)){tmp[1:(op-1)]=otmp[1:(op-1)] }
if(op!=1 && (op+1)!=length(tmp)){
tmp[(op+2):length(tmp)]=otmp[(op+1):(length(tmp)-1)]
tmp[1:(op-1)]=otmp[1:(op-1)]}
rownames(mtx)<-c(tmp)
colnames(mtx)<-c(tmp)
mtx<- fillspcmtx()
otmp<-tmp
omtx<-mtx
#above generate sequence and cov. matrix for speciation.
}else{# op = 0 means that we have detected the hybridize event.
ins<-(bmtp[q])#grab the insertion element, ins will be used in the fillhydmtx function.
insp<-ins_fcn(ins,otmp)#catch the position for insertion.
tmp[insp]<-ins #insert the hyd element.
tmp[(insp+1):length(tmp)]=otmp[insp:(length(tmp)-1)]
tmp[1:(insp-1)]=otmp[1:(insp-1)]
rownames(mtx)<-c(tmp)
colnames(mtx)<-c(tmp)
diamtx<-var_fcn()
mtx<- fillhydmtx()
mtx<-mtx+t(mtx)+diamtx
#fill in the diagonal elements.
otmp<-tmp
omtx<-mtx
#above generate the sequnce and fill in the value into matrix for hybrdization.
ckins<-TRUE #since we did an insertion, the next step is to replace 3 elements.
}#end of the length(op)!=0 if-else.
if(ckins){#replace 3 elements in tmp sequence.
tmp<-array(0,c(length(tmp)))
tmp[which(otmp==ins)]<- rept[loop-1,1] # replaced with hyd element.
tmp[which(otmp == bmtp[loop])] = rept[loop,which(rept[loop,]!=ins)]
tmp[which(otmp == bmtp[loop+1])] = rept[loop+1,which(rept[loop+1,]!=ins)]
#replace 3 new nodes.
tx1<-which(otmp==ins)
tx2<-which(otmp == bmtp[loop])
tx3<-which(otmp == bmtp[loop+1])
for(i in 1:length(tmp)){
if (i != tx1 && i!=tx2 && i!=tx3){
tmp[i]=otmp[i]
}
}
otmp<-tmp
rownames(mtx)<-c(tmp)
colnames(mtx)<-c(tmp)
mtx<-mtx+diag(branchlength[length(tmp)-1],c(length(tmp)) )
omtx<-mtx
ckins<-FALSE
loop<-loop+2
}#end of replace 3 elements
}#end of while loop
if(sum(tipnames%in%tmp)!=nleaves){#catches the last speciation event.
tmp<-tipnames
mtx<-matrix(0,nrow=length(tmp),ncol=length(tmp))
rownames(mtx)<-c(tmp)
colnames(mtx)<-c(tmp)
mtx<-fillspcmtx()
}#end of if (sum(tipnames%in%tmp)!=nleaves).
mtx<-mtx*sigma_sq
hybrid.Index<-hybrid.node(ancdes.array,nleaves) #extra burst for hybrid
for(i in hybrid.Index){
mtx[i,i]<-mtx[i,i]+sigma.H_sq
}
if(model.Index==1){rm(bt)}
if(model.Index==2){rm(sigma.H_sq)}
if(model.Index==3){rm(bt);rm(sigma.H_sq)}
return(mtx)
}#end of cov_mtx
hybrid.node<-function(ancdes.array,nleaves){
hyd.sigma_h<-c()
for (i in 1:dim(ancdes.array)[1]){
if(ancdes.array[i,2]==ancdes.array[i,3]){
hyd.idx<-as.numeric(unlist(strsplit(ancdes.array[i,2],"X"))[2])
if ( hyd.idx<=nleaves){
hyd.sigma_h<-c(hyd.sigma_h,hyd.idx )}else{
hyd.speciation<-which(ancdes.array[,1]==ancdes.array[i,2])
candidate.hyd.des<-as.numeric(unlist(strsplit(ancdes.array[hyd.speciation,],"X")[2:3]))[c(2,4)]
for(hyd.des in candidate.hyd.des){
if(hyd.des<=nleaves){
hyd.sigma_h<-c(hyd.sigma_h,hyd.des )}}
}
}
}
return(hyd.sigma_h )
}
#sunflower 11 taxa (3 hybrids) Gross and Riesberg paper
x.tre<-c("(((1:17,2:17)21:17,(3:19,(((4:6,(5:6)13:0)12:5,(6:11)16:0)15:4,(7:15)19:0)18:4)22:15)26:6,(((19:0,(16:0,(13:0,8:6)14:5)17:4)20:5,9:20)23:6,(10:23,11:23) 24:3)25:14)27:0;")
res<-newick2phylog(x.tre)
ancdes.array<-getntn(res)
branchlength<-getbrnlen(res)
branchlength<-branchlength[-length(branchlength)]
tipnames<-sort(names(res$droot[which(res$droot==max(res$droot))]))
nleaves<-length(tipnames)
n<-nleaves
model.Index<-2
load("regsunflower.RData")check for positive definiteness of Gtau
The script in question is designed to test whether the variance-covariance matrix for the network model \(G_\hat{\tau}\) is positive definite. This test is performed after the MLE estimate is given, and it is important to ensure this regularity condition is met, as it guarantees that the regression parameters are normally distributed. This, in turn, allows for maximum likelihood analysis to be performed.
The results of the script can be found at the bottom of the webpage, and they indicate that the positive definite property is satisfied for all 12 data sets used in this analysis.
check point.
To ensure the likelihood satisfies the regularity condition, it is necessary to check whether \(G_{\hat{\tau}}\) is positive definite. Please see the following
h<-0.5
sigma_sq<-1
sigma.H_sq<-0
for(btIndex in 1:dim(MLE.array)[1]){
tau<-MLE.array[btIndex,1]
x.mle=c(tau,h,sigma_sq,sigma.H_sq)
V_hat<-cov_mtx(x=x.mle,branchlength=branchlength,ancdes.array=ancdes.array,nleaves=nleaves,tipnames=tipnames,model.Index=model.Index)
print(paste("data", btIndex,". Is positve definite ?",sep=""))
print(is.positive.definite(V_hat, tol=1e-8))
}[1] "data1. Is positve definite ?"
[1] TRUE
[1] "data2. Is positve definite ?"
[1] TRUE
[1] "data3. Is positve definite ?"
[1] TRUE
[1] "data4. Is positve definite ?"
[1] TRUE
[1] "data5. Is positve definite ?"
[1] TRUE
[1] "data6. Is positve definite ?"
[1] TRUE
[1] "data7. Is positve definite ?"
[1] TRUE
[1] "data8. Is positve definite ?"
[1] TRUE
[1] "data9. Is positve definite ?"
[1] TRUE
[1] "data10. Is positve definite ?"
[1] TRUE
[1] "data11. Is positve definite ?"
[1] TRUE
[1] "data12. Is positve definite ?"
[1] TRUE