rm(list=ls())
setwd("~/Dropbox/NSCT/112年計畫113年八月/phygarch11rate/rcode/simulation/")
library(ape)
library(phytools)
library(geiger)
library(TreeSim)
library(tseries)
library(phyclust)
#source("~/Dropbox/MOST_NSC/111年計畫112年八月/risingPp/fcliu/rcode/")
taxasize.array <- c(16, 32, 64, 128) #c(22)
treetype.array <- c("balanced")#,"left","coalescent","bdtree") #c("bdtree")
root.y <- 0.2
root.sig <- 0.0001
root.eps <- 0
omega.true <- 1.5e-1
alpha.true <- 0.00852
beta.true <- 1e-7
Theta.true<-c(omega.true,alpha.true,beta.true)
names(Theta.true)<-c("omega","alpha","beta")
for(taxasizeIndex in 1:length(taxasize.array)){
# taxasizeIndex<-1
taxasize<-taxasize.array[taxasizeIndex]
taxasize
for(treetypeIndex in 1:length(treetype.array)){
# treetypeIndex<-1
treetype<-treetype.array[treetypeIndex]
treetype
if(treetype=="balanced"){tree<-compute.brlen(stree(taxasize,type="balanced"))}
if(treetype=="left"){tree<-compute.brlen(stree(taxasize,type="left"))}
if(treetype=="coalescent"){tree<-rcoal(taxasize)}
if(treetype=="bdtree"){
brith_lambda=runif(1,0.01,0.1)
death_mu=runif(1,0,brith_lambda)
tree<-sim.bd.taxa(n=taxasize,numbsim=1,lambda=brith_lambda,mu=death_mu,complete = FALSE, stochsampling = TRUE)[[1]]
}
tree= reorder(tree,"postorder")
# plot(tree)
#vcv(tree)
anc =tree$edge[,1]
des =tree$edge[,2]
tree$edge.length<-tree$edge.length/get.rooted.tree.height(tree,tol = .Machine$double.eps^0.5)
treelength<-tree$edge.length
treelength
N=dim(tree$edge)[1]
ynodestate= array(NA,c(2*taxasize-1))
signodestate= array(NA,c(2*taxasize-1))
epsnodestate= array(NA,c(2*taxasize-1))
ynodestate[taxasize+1]<-root.y
signodestate[taxasize+1]<-root.sig
epsnodestate[taxasize+1]<-root.eps
# need to use true params to simulate trait first
# This simulate create crazy big number of sims value
for(Index in N:1){
# Index<-N
signodestate[des[Index]] = Theta.true['omega'] + Theta.true['alpha']*(epsnodestate[anc[Index]])^2 + Theta.true['beta']*(signodestate[anc[Index]])^2
epsnodestate[des[Index]]<-signodestate[des[Index]]*rnorm(1,0,sd=sqrt(treelength[Index]))
ynodestate[des[Index]]=ynodestate[anc[Index]]+epsnodestate[des[Index]]
}# end of tree traversal
signodestate # small variation for the sigma shall be set up for sims good
epsnodestate # the epsilon is too. that would contributed to the next trait node values
ynodestate # we may suppose that the node value shall not vary tremendous
#plot(signodestate)
#summary(signodestate)
#mean(signodestate)
#mean(ynodestate)
# gest<-garch(ynodestate, order = c(1, 1))
# Theta<-gest$coef
# names(Theta)<-c("omega","alpha","beta")
# Theta.garch<-Theta
# Theta.garch
# hyper_omega <- 3*Theta["omega"]
# hyper_alpha <- 3*Theta["alpha"]
# hyper_beta <- 3*Theta["beta"]
# here now we have sim data from the true parameters, we need to access the performance of model.
# use the tip for inference ?
trait<-ynodestate[1:taxasize]
names(trait)<-tree$tip.label
# fitbm = fitContinuous(phy=tree,dat=trait)
# sig0 = sqrt(fitbm$opt$sigsq)
# y0 = fitbm$opt$z0
#first compute ell0
names(trait)=tree$tip.label
# fitbm=fitContinuous(phy=tree,dat=trait)
# sig0=sqrt(fitbm$opt$sigsq)
# y0=fitbm$opt$z0
# eps0<-(trait-y0)
#ell0 <- -fitbm$opt$lnL
#ell0 <- -1/2*log(2*pi)-1/2*log(root.sig)-1/2/root.sig*sum(root.eps^2)
#ell0<- 0# -Inf # can we set up inf? I think so, because it would make the ell-ell0 to Inf,
# then the min between 1 and the ratio of the likelihood is 1, and the first proposal is accepted always,
# the, you replace the new likelihood ell0 with the ell shall be ok
eps0<-(trait-root.y)
ell0 <- -1/2*log(2*pi)-1/2*log(root.sig)-1/2/root.sig*sum(root.eps^2)
for(Index in N:1){
# if(signodestate[des[Index]]<0){signodestate[des[Index]]<-1}
ell0 = ell0 + 1/2*(-log(2*pi)-log(signodestate[des[Index]])-epsnodestate[des[Index]]/signodestate[des[Index]])
# print(ell)
# print(c(Index,signodestate[des[Index]]))
}
ell0
#################
#################
n_iter<-1e6
n_burnin<-n_iter/2
ell.array<-array(NA,n_iter)
adapt=FALSE #RAM algorithm
accept<-array(0,n_iter)
Theta.vec<-array(NA,c(n_iter,length(Theta.true)))
#install.packages("tseries")
Theta.vec[1,]<-Theta.true
colnames(Theta.vec)<-c("omega","alpha","beta")
accept[1]<-1
ell0
# Smtx<-diag(length(Theta))
for (simIndex in 2:n_iter){
if(simIndex %% 5000 == 0){print(
paste(simIndex," out of ", n_iter, "interation",sum(accept), " accepted, acc rate = ", sum(accept)/n_iter,sep=""))}
# simIndex<- 2
# omega <- runif(1,0,hyper_omega) #0.00005
# alpha <- runif(1,0,hyper_alpha) #0.085
# beta <- runif(1,0,hyper_beta) #0.01
### Here garch estimate
# omega<-rnorm(1,mean=Theta.garch["omega"],sd=Theta.garch["omega"]/4)
# alpha<-rnorm(1,mean=Theta.garch["alpha"],sd=Theta.garch["alpha"]/4)
# beta<-rnorm(1,mean=Theta.garch["beta"],sd=Theta.garch["beta"]/4)
Theta.current<-Theta.vec[simIndex-1,]
### Here true estimate
# omega.prop<-Theta.current["omega"]+rnorm(1,mean=0,sd=Theta.true["omega"]/4)
### Here true estimate
omega.prop<-omega.true#Theta.current["omega"]+rnorm(1,mean=0,sd=Theta.true["omega"]/4)
alpha.prop<-abs(Theta.current["alpha"]+rnorm(1,mean=0,sd=Theta.true["alpha"]/2))
beta.prop<- abs(Theta.current["beta"]+rnorm(1,mean=0,sd=Theta.true["beta"]/2))
Theta.prop <- c(omega.prop,alpha.prop,beta.prop)
names(Theta.prop)<-c("omega","alpha","beta")
Theta.prop
Theta.current
# Uvec<-matrix(rnorm(length(Theta),mean=Theta,sd=Theta/2),ncol=1)
# Theta_new<-c(Theta+Smtx%*%Uvec)
# names(Theta_new)<-c("omega","alpha","beta")
# Theta_new
###########################
## simulate nodes values ##
###########################
ynodestate= array(NA,c(2*taxasize-1))
signodestate= array(NA,c(2*taxasize-1))
epsnodestate= array(NA,c(2*taxasize-1))
ynodestate[taxasize+1]<-root.y#y0
signodestate[taxasize+1]<-root.sig#sig0
epsnodestate[taxasize+1]<-root.eps#mean(trait-y0)
for(Index in N:1){
# Index<-N
signodestate[des[Index]]=Theta.prop['omega']+Theta.prop['alpha']*(epsnodestate[anc[Index]])^2 + Theta.prop['beta']*(signodestate[anc[Index]])^2
epsnodestate[des[Index]]<-signodestate[des[Index]]*rnorm(1,0,sd=sqrt(treelength[Index]))
ynodestate[des[Index]]=ynodestate[anc[Index]]+epsnodestate[des[Index]]
}# end of tree traversal
signodestate
epsnodestate
ynodestate
########################
## compute likelihood ##
########################
ynodestate[1:taxasize]=trait
eps0<-(trait-root.y)
ell <- -1/2*log(2*pi)-1/2*log(root.sig)-1/2/root.sig*sum(root.eps^2)
for(Index in N:1){
# if(signodestate[des[Index]]<0){signodestate[des[Index]]<-1}
ell = ell + 1/2*(-log(2*pi)-log(abs(signodestate[des[Index]]))-epsnodestate[des[Index]]/signodestate[des[Index]])
# print(ell)
# print(c(Index,signodestate[des[Index]]))
}
ell
ell0
ell.array[simIndex]<-ell
delta<- min(1,exp(ell-ell0)) # here the new proposal generate a new likelihood, the ratio greater 1, it just take 1(so the minimum) between the two. if less than 1, then let a uniform random variable to determine the acceptance by the comparison. the chance is there for being accepted aor rejected.
delta
if( runif(1)< delta){ # here is the acceptance . the the sample is updated
# print(paste(simIndex," accept ", "ell=",round(ell,2),", ell0=",round(ell0,2),sep=""))
accept[simIndex]<-1
#accept[1:2]
Theta.vec[simIndex,]<- Theta.prop
Theta.current<-Theta.prop
ell0<-ell
}else{
# print("reject proposal")
# print(paste(simIndex," reject ", "ell=",round(ell,2),", ell0=",round(ell0,2),sep=""))
Theta.vec[simIndex,]<- Theta.vec[simIndex-1,]
accept[simIndex]<-0
}
head(Theta.vec)
if(adapt==TRUE & simIndex <= n_burnin){
S=ramcmc::adapt_S(S=Smtx,u=Uvec,delta,simIndex-1)
}
}#end of n_iter
# par(mfrow=c(1,3))
# plot(Theta.vec[1:n_iter,"omega"],main="omega")
# plot(Theta.vec[1:n_iter,"alpha"],main="alpha")
# plot(Theta.vec[1:n_iter,"beta"],main="beta")
#
# output<-list(Theta = Theta.vec[(n_burnin + 1):n_iter, ], acceptance_rate = sum(accept[(n_burnin + 1):n_iter]) / (n_iter - n_burnin))
# output$acceptance_rate
#
# plot(output$Theta[,"omega"],main="omega")
# plot(output$Theta[,"alpha"],main="alpha")
# plot(output$Theta[,"beta"],main="beta")
#
# head(output)
# apply(output$Theta,2,mean,na.rm=T)
# apply(output$Theta,2,median,na.rm=T)
# apply(output$Theta,2,sd,na.rm=T)
#
# dim(output$Theta)
# apply(output$Theta,2,summary)
# save.image(paste(treetype,"taxasize", taxasize, ".rda",sep=""))
#Theta
#if(FALSE){
#
# plot(output$Theta[,"omega"])
# plot(output$Theta[,"alpha"])
# plot(output$Theta[,"beta"])
#
# par(mfrow=c(1,3))
# hist(output$Theta[,"omega"])
# hist(output$Theta[,"alpha"])
# hist(output$Theta[,"beta"])
# #}
# acf(output$Theta)
# plot(density(output$Theta[,"omega"]))
# plot(density(output$Theta[,"alpha"]))
# plot(density(output$Theta[,"beta"]))
#save.image(paste(treetype,"taxa",taxasize,".rda",sep=""))
}#treetypeIndex
}#taxasizeIndexModel Script with Simulation
Various type of trees: star tree, balanced tree, left tree, coalescent tree, birth-death tree will be used for assessment using taxa \(8, 16, 32, 64, 128, 256\) and \(512\).
Two sets of priors will be used: the informative prior and the flat prior.
For each type of tree, for each data and for each prior set, a thousand replicates of the sample will be generated to assess the posterior sample distribution.
load("~/Dropbox/NSCT/112年計畫113年八月/phygarch11rate/rcode/simulation/balancedtaxasize128.rda")
par(mfrow=c(1,3))
plot(Theta.vec[1:n_iter,"omega"],main="omega")
plot(Theta.vec[1:n_iter,"alpha"],main="alpha")
plot(Theta.vec[1:n_iter,"beta"],main="beta")
output<-list(Theta = Theta.vec[(n_burnin + 1):n_iter, ], acceptance_rate = sum(accept[(n_burnin + 1):n_iter]) / (n_iter - n_burnin))
output$acceptance_rate[1] 0.3086plot(output$Theta[,"omega"],main="omega")
plot(output$Theta[,"alpha"],main="alpha")
plot(output$Theta[,"beta"],main="beta")
#head(output)
apply(output$Theta,2,mean,na.rm=T) omega alpha beta
1.500000e-01 2.243686e-01 1.959763e-06 apply(output$Theta,2,median,na.rm=T) omega alpha beta
1.500000e-01 2.279207e-01 1.748528e-06 apply(output$Theta,2,sd,na.rm=T) omega alpha beta
0.000000e+00 2.862598e-02 6.436357e-07 dim(output$Theta)[1] 5000 3apply(output$Theta,2,summary) omega alpha beta
Min. 0.15 0.1609715 9.299612e-07
1st Qu. 0.15 0.1992346 1.475306e-06
Median 0.15 0.2279207 1.748528e-06
Mean 0.15 0.2243686 1.959763e-06
3rd Qu. 0.15 0.2468118 2.396948e-06
Max. 0.15 0.2920027 3.709344e-06#Theta
#if(FALSE){
#
par(mfrow=c(1,3))
plot(output$Theta[,"omega"],type="l")
plot(output$Theta[,"alpha"],type="l")
plot(output$Theta[,"beta"],type="l")
#
par(mfrow=c(1,3))
hist(output$Theta[,"omega"])
hist(output$Theta[,"alpha"])
hist(output$Theta[,"beta"])