---
title: "check for positive definiteness of Gtau"
format: html
editor: visual
---

The following script check the positive definiteness of 
the variance covariance matrix for the network model 
$G_\hat{\tau}$ after given the MLE estimate in order to ensure to regularity condition that guarantee the normality condition for the regression parameters(hence the inference can be proceed through maximum likelihood analysis).

Results (see the the bottom of this webpage page) shows the positive definite property can be met for all the 12  data set used here. 

```{r}
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 the $G_{\hat{\tau}}$ is positive definite to ensure the regularity condition for the likelihood. 

```{r}
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))
}
```





