rm(list=ls())
library(ape)
library(phytools)
## Loading required package: maps
tree<-"((z:1,(y:0.42,x:0.42):0.58):0);"
phy<-read.newick(text=tree)

tree2<-"((y:1,(z:0.72,x:0.72):0.28):0);"
phy2<-read.newick(text=tree2)

phy$tip.label<-c("   y","   x","   z")

phy2$tip.label<-c("   z","   y","   z")

colors1<-c("red","orange","blue","black")

colors2<-c("blue","orange","red","black")

phy$edge.length<-phy$edge.length*100
phy2$edge.length<-phy2$edge.length*100

op<-par(mfrow=c(1,2),
        oma=c(2,2,0,0)+0.1,
        mar=c(1,1,1,1)+0.1        
)

plot(phy,edge.width=5,cex=2,edge.color=rev(colors1))
#tiplabels(pch=21, col="black", adj=1, bg="black", cex=2)
axisPhylo(1, root.time = NULL, backward = F)

plot(phy2,edge.width=5,cex=2,edge.color=rev(colors2))
#tiplabels(pch=21, col="black", adj=1, bg="black", cex=2)
axisPhylo(1, root.time = NULL, backward = F)

ones<-matrix(c(1,1,1),ncol=1)
X1<-matrix(c(1.2,1.5,0.5),ncol=1)
C1<-vcv(phy)
mu1hat<- t(ones)%*%solve(C1)%*%X1/sum(solve(C1))
sigsq1hat <- t(X1-c(mu1hat)*ones)%*%solve(C1)%*%(X1-c(mu1hat)*ones)/3

negLike1<-3/2*log(2*pi)+0.5*log(det(c(sigsq1hat)*C1)) + 1/2/sigsq1hat*t(X1-c(mu1hat)*ones)%*%solve(C1)%*%(X1-c(mu1hat)*ones)
negLike1
##          [,1]
## [1,] 2.693285
X2<-c(1.2,0.5,1.5)
C2<-vcv(phy2)
mu2hat<- t(ones)%*%solve(C2)%*%X2/sum(solve(C2))
sigsq2hat <- t(X2-c(mu2hat)*ones)%*%solve(C2)%*%(X2-c(mu2hat)*ones)/3
negLike2<-3/2*log(2*pi)+0.5*log(det(c(sigsq2hat)*C2))+1/2/sigsq2hat*t(X2-c(mu2hat)*ones)%*%solve(C2)%*%(X2-c(mu1hat)*ones)
negLike2
##          [,1]
## [1,] 2.072889
c(mu1hat,sigsq1hat)
## [1] 1.088268156 0.004042742
c(mu2hat,sigsq2hat)
## [1] 1.078048780 0.002396116