---
title: "Model Script with Simulation"
format: html
editor: visual
---

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.

```{r,eval=FALSE}
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
}#taxasizeIndex

```


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

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)
#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"])
```