diff --git a/brnlenMCMC.py b/brnlenMCMC.py index 05ebc39..26a5d96 100644 --- a/brnlenMCMC.py +++ b/brnlenMCMC.py @@ -169,7 +169,7 @@ def replaedglen(): def mcmcbrn(postorder, patterns, rates): nodes = readnewick(treenewick()) - run = 5 + run = 2 mcmc = 0 output = os.path.join('brnlenMCMC_results.txt') newf = open(output, 'w') @@ -194,34 +194,32 @@ def mcmcbrn(postorder, patterns, rates): print # print i, nd.number - print 'current', - for j in nodes: - print j.edgelen, +# print 'current', + print 'current_node.edgelen=', nodes[i].number, nodes[i].edgelen # print preedgelen = nodes[i].edgelen # print 'nodes[i]=', nodes[i].edgelen # prelike = prepareTree(postorder, patterns, rates) - prelike = prepareTree(nodes, patterns, rates) - preprior = (-nodes[i].edgelen/b)-(log(b)) - pre_ln_posterior = prelike + preprior + currentlike = prepareTree(nodes, patterns, rates) + currentprior = (-nodes[i].edgelen/b)-(log(b)) + current_ln_posterior = currentlike + currentprior u = random.random() m = exp(0.2*(u-0.5)) nodes[i].edgelen = preedgelen*m + print 'proposed_node.edgelen=',nodes[i].number, nodes[i].edgelen # print 'new_edgelen=', nodes[i].edgelen - postprior = (-nodes[i].edgelen/b)-(log(b)) - postlike = prepareTree(nodes, patterns, rates) - post_ln_posterior = postlike + preprior + proposedprior = (-nodes[i].edgelen/b)-(log(b)) + proposedtlike = prepareTree(nodes, patterns, rates) + proposed_ln_posterior = proposedtlike + proposedprior print -# print 'preprior=', preprior - print 'post_ln_posterior=', post_ln_posterior - print 'pre_ln_posterior=', pre_ln_posterior -# print 'postlike=', postlike + print 'current prior, currentlike=', currentprior, currentlike + print 'proposed prior, proposedlike=', proposedprior, proposedtlike print hastings_ratio = log(m) print 'hastings_ratio', hastings_ratio - logR = post_ln_posterior - pre_ln_posterior + hastings_ratio + logR = proposed_ln_posterior - current_ln_posterior + hastings_ratio print 'logR=', logR u2 = random.random() print 'log-u2=', log(u2) @@ -231,11 +229,7 @@ def mcmcbrn(postorder, patterns, rates): else: nodes[i].edgelen = preedgelen print 'failed accept new proposal..' - # print 'nodes[i]-b=============================', nodes[i].edgelen - # for i in nodes: - # print 'new', i.edgelen - - # print 'pre, post=', preprior, postprior + newf.write('%s\t'%(mcmc+1)) print print 'proposed=', @@ -250,28 +244,6 @@ def mcmcbrn(postorder, patterns, rates): print '------------------------------------------------------------------------------' -# b = 0.1 -# preedgelen = node.edgelen - - -# preprior = (-(node.edgelen)/b)-(log(b)) -# pre_ln_posterior = preloglike + log(preprior) -# -# u = random.random() -# m = exp(0.2*(u-0.5)) -# newedgelen = node.edgelen*m -# postprior = (-(node.newedgelen)/b)-(log(b)) -# post_ln_posterior = postloglike + log(postprior) -# -# hastings_ratio = log(m) -# logR = post_ln_posterior - pre_ln_posterior + hastings_ratio -# -# u2 = random.random() -# if log(u2) < logR: -# node.edgelen = newedgelen -# else: -# node.edgelen = node.edgelen -# def treenewick(): @@ -490,10 +462,10 @@ def calcExpectedHeight(num_species, mu_over_s): if __name__ == '__main__': -# random_seed = 146549 # 7632557, 12345 + random_seed = 146549 # 7632557, 12345 number_of_species = 5 mutation_speciation_rate_ratio = 0.4 # 0.689655172 # yields tree height 1 for 6 species -# random.seed(random_seed) + random.seed(random_seed) species_tree_root = yuleTree(number_of_species, mutation_speciation_rate_ratio) # print '#########' # print species_tree_root diff --git a/treeLike_JCGAMMA.py b/treeLike_JCGAMMA.py index c14e0ce..4011773 100644 --- a/treeLike_JCGAMMA.py +++ b/treeLike_JCGAMMA.py @@ -160,7 +160,7 @@ def allocatePartial(node, patterns, rates): def treenewick(): script_dir = os.path.dirname(os.path.realpath(sys.argv[0])) - path = os.path.join(script_dir, 'tree.tre') + path = os.path.join(script_dir, 'tree2.tre') with open(path, 'r') as content: newick = content.read() return newick