-
Notifications
You must be signed in to change notification settings - Fork 0
/
Traditional_MH.R
165 lines (136 loc) · 6.88 KB
/
Traditional_MH.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
#### libraries loaded ####
library(ape)
library(phylobase)
library(mvtnorm)
library(MASS)
#### hyperParam ####
source("./Downloads/MH_Func.R")
args <- commandArgs(TRUE)
print(args)
c_rate=c_shape=2
s_rate=s_shape=1
c_vec<-as.integer(c(seq(0.3,1,0.1),seq(1.5,3,0.5))*10)/10
sigma_vec<-seq(0.5,2,0.5)
obsData_idx<-as.numeric(args[1])
true_c<-as.integer(as.numeric(args[2])*10)/10
true_sigma2<-as.integer(as.numeric(args[3])*10)/10
chain_idx<-as.numeric(args[4])
c_idx<-which(c_vec==true_c)
s_idx<-which(sigma_vec==true_sigma2)
N=100
d=10
setwd('~/Research_Veera/PYDT/MCMC_result/DDT_MH')
#### Obs Data ####
set.seed(202010023+ceiling(obsData_idx*100000)+ceiling(chain_idx*10000)+ceiling(c_idx*100)+ceiling(s_idx*100))
obsPhylo4d_lt<-readRDS(paste0("~/Research_Veera/PYDT/ABC_simData/DDT/N100d10/obsData/obs_phylo4d_lt",obsData_idx,".RDS"))
#obsPhylo4d_lt<-readRDS(paste0("~/Research_Veera/PYDT/ABC_simData/DDT/N100d10/obsData/obs_phylo4d_lt1.RDS"))
obsDf<-obsPhylo4d_lt[[(c_idx-1)*length(sigma_vec)+s_idx]]@data[1:N,1:d]
rownames(obsDf)<-paste0("N",1:N)
if(as.integer(obsPhylo4d_lt[[(c_idx-1)*length(sigma_vec)+s_idx]]@data[1,"c"]*10)/10!=true_c |
as.integer(obsPhylo4d_lt[[(c_idx-1)*length(sigma_vec)+s_idx]]@data[1,"sigma2"]*10)/10!=true_sigma2){
stop("different true param with df")
}
#### Init ####
options(warn=-1)
c_old<-rgamma(1,shape = c_shape, rate = c_rate)
s2_old<-1/rgamma(1,shape = s_shape, rate = s_rate)
method_vec=c("ward.D", "ward.D2", "single", "complete", "mcquitty")
obs_hclu_ape<-as.phylo(hclust(dist(obsDf),method=method_vec[chain_idx]))
root_edge0<-div_time(1,0,c_old,0,0)
obs_hclu_ape<-add_root(tr=obs_hclu_ape, root_edge = root_edge0, rootLabel = paste0("N",N+1), tipLabel = paste0("N","N+2"))
obs_hclu_phylo4<-phylo4(obs_hclu_ape)
obs_hclu_phylo4_totalLen<-nodeHeight(obs_hclu_phylo4,"N1","root")
edgeLength(obs_hclu_phylo4)<-edgeLength(obs_hclu_phylo4)/obs_hclu_phylo4_totalLen
obs_hclu_phylo4d<-phylo4d(obs_hclu_phylo4,tip.data=obsDf)
nodeHeight(obs_hclu_phylo4d,"N1","root")
nodeLabels(obs_hclu_phylo4d)<-paste0("N",(N+1):(2*N))
saveRDS(obs_hclu_phylo4d,paste0("./MCMC_Compare/N100d10_Replicate/DefaultMH/init_phylo4d_lt",obsData_idx,"_c",true_c,"_s",true_sigma2,"_ch",chain_idx,".RDS"))
#### MH Alg ####
#MCMC_num<-30000
#MCMC_num<-10000
MCMC_num<-5000
#MCMC_num<-3
print(paste0("total MCMC number: ",MCMC_num))
llh_chk<-data.frame(idx=numeric(0),
Accept=numeric(0),
llh_Prop=numeric(0),
q_ToOld=numeric(0),
llh_old=numeric(0),
q_ToNew=numeric(0),
c=numeric(0),
sigma2=numeric(0))
tree_lt<-list()
MRCA_mat_lt<-list()
trStruct_lt<-list()
options(warn=-1)
system.time({
for(i in 1:MCMC_num){
print(i)
if(i==1){
old_phylo4d<-obs_hclu_phylo4d
}
old_phylo4<-extractTree(old_phylo4d)
rnd_Detach<-Random_Detach(tr_phylo4 = old_phylo4)
prop_newTree<-attach_subTree(subTree = rnd_Detach$subTr, rmnTree = rnd_Detach$rmnTr,
old_divTime = rnd_Detach$divT, old_subTr_paLab = rnd_Detach$paNodeLab,
c_ = c_old, theta_=0, alpha_=0)
while(nTips(prop_newTree$new_phylo4)!=N){
prop_newTree<-attach_subTree(subTree = rnd_Detach$subTr, rmnTree = rnd_Detach$rmnTr,
old_divTime = rnd_Detach$divT, old_subTr_paLab = rnd_Detach$paNodeLab,
c_ = c_old, theta_=0, alpha_=0)
}
q_prob<-prop_log_prob(orgTree_phylo4 = old_phylo4, rmnTree = rnd_Detach$rmnTr,
old_detach_PaTime = rnd_Detach$paDivT, old_detach_PaLab = rnd_Detach$paNodeLab, old_detach_divLab = rnd_Detach$divLab,
new_divTime = prop_newTree$new_divT, new_AttachRoot = prop_newTree$attach_root, new_AttachTo = prop_newTree$attach_to,
c_=c_old)
new_phylo4d<-phylo4d(prop_newTree$new_phylo4,tip.data=obsDf)
if(i==1){
oldTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = old_phylo4d)
newTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = new_phylo4d)
}else{
oldTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = old_phylo4d,
trStruct_old = trStruct_lt[[(i-1)]] , tre_sigma_old = MRCA_mat_lt[[(i-1)]])
newTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = new_phylo4d)
}
while(length(newTree_info)==1 | nTips(prop_newTree$new_phylo4)!=N){
prop_newTree<-attach_subTree(subTree = rnd_Detach$subTr, rmnTree = rnd_Detach$rmnTr,
old_divTime = rnd_Detach$divT, old_subTr_paLab = rnd_Detach$paNodeLab,
c_ = c_old, theta_=0, alpha_=0)
q_prob<-prop_log_prob(orgTree_phylo4 = old_phylo4, rmnTree = rnd_Detach$rmnTr,
old_detach_PaTime = rnd_Detach$paDivT, old_detach_PaLab = rnd_Detach$paNodeLab, old_detach_divLab = rnd_Detach$divLab,
new_divTime = prop_newTree$new_divT, new_AttachRoot = prop_newTree$attach_root, new_AttachTo = prop_newTree$attach_to,
c_=c_old)
new_phylo4d<-phylo4d(prop_newTree$new_phylo4,tip.data=obsDf)
if(i==1){
oldTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = old_phylo4d)
newTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = new_phylo4d)
}else{
oldTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = old_phylo4d,
trStruct_old = trStruct_lt[[(i-1)]] , tre_sigma_old = MRCA_mat_lt[[(i-1)]])
newTree_info<-DDT_log_llh(c_=c_old, cmn_sigma2 = s2_old, tr_phylo4d = new_phylo4d)
}
}
llh_chk[i,-2]<-c(i,newTree_info$res_log_llh,q_prob["q_RmnToOld"],
oldTree_info$res_log_llh,q_prob["q_RmnToNew"],
c_old, s2_old)
u<-runif(1)
alpha_Accp<-llh_chk[i,"llh_Prop"]+llh_chk[i,"q_ToOld"]-llh_chk[i,"llh_old"]-llh_chk[i,"q_ToNew"]
if(u<exp(alpha_Accp)){
llh_chk[i,2]<-1
tree_lt[[i]]<-new_phylo4d
MRCA_mat_lt[[i]]<-newTree_info$tre_sigma
trStruct_lt[[i]]<-newTree_info$trStruct
old_phylo4d<-new_phylo4d
}else{
llh_chk[i,2]<-0
tree_lt[[i]]<-old_phylo4d
MRCA_mat_lt[[i]]<-oldTree_info$tre_sigma
trStruct_lt[[i]]<-oldTree_info$trStruct
}
c_old<-postSamp_c(shape_prior = c_shape, rate_prior = c_rate, trStruct_old = trStruct_lt[[i]])
s2_old<-1/postSamp_s(shape_prior = s_shape, rate_prior = s_rate, MRCA_mat_old = MRCA_mat_lt[[i]])
}
})
saveRDS(llh_chk,paste0("./MCMC_Compare/N100d10_Replicate/DefaultMH/llh_lt",obsData_idx,"_c",true_c,"_s",true_sigma2,"_ch",chain_idx,".RDS"))
saveRDS(tree_lt,paste0("./MCMC_Compare/N100d10_Replicate/DefaultMH/tree_lt",obsData_idx,"_c",true_c,"_s",true_sigma2,"_ch",chain_idx,".RDS"))
options(warn=0)