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