Cleanup clocks

tutorials/clocks/index.md

@@ -342,25 +342,29 @@ still depend on variable initialized in different files.
 
 The clock-rate parameter is a stochastic node from a gamma distribution.
 
+{% snippet scripts/m_GMC_bears.Rev %}
     clock_rate ~ dnGamma(2.0,4.0)
     moves.append( mvScale(clock_rate,lambda=0.5,tune=true,weight=5.0) )
+{% endsnippet %}
 
 ***The Sequence Model and Phylogenetic CTMC***
 
 Specify the parameters of the GTR model and the moves to operate on
 them.
-```
+
+{% snippet scripts/m_GTR.Rev %}
     sf ~ dnDirichlet(v(1,1,1,1))
     er ~ dnDirichlet(v(1,1,1,1,1,1))
     Q := fnGTR(er,sf)
     moves.append( mvSimplexElementScale(er, alpha=10.0, tune=true, weight=3.0) )
     moves.append( mvSimplexElementScale(sf, alpha=10.0, tune=true, weight=3.0) )
-```
+{% endsnippet %}
+
 And instantiate the phyloCTMC.
-```
+{% snippet scripts/m_GMC_bears.Rev %}
     phySeq ~ dnPhyloCTMC(tree=timetree, Q=Q, branchRates=clock_rate, nSites=n_sites, type="DNA")
     phySeq.clamp(D)
-```
+{% endsnippet %}
 This is all we will include in the global molecular clock model file.
 
 Save and close the file called in the `scripts` directory.
@@ -378,30 +382,30 @@ the `scripts` directory.
 
 *Load Sequence Alignment* — Read in the sequences and initialize
 important variables.
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
     D <- readDiscreteCharacterData(file="data/bears_irbp.nex")
     n_sites <- D.nchar()
     mi = 1
-
+{% endsnippet %}
 *The Calibrated Time-Tree Model* — Load the calibrated tree model from
 file using the `source()` function. Note that this file does not have
 moves that operate on the tree topology, which is helpful when you plan
 to estimate the marginal likelihoods and compare different relaxed clock
 models.
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
     source("scripts/m_BDP_bears.Rev")
-
+{% endsnippet %}
 *Load the GMC Model File* — Source the file containing all of the
 parameters of the global molecular clock model. This file is called .
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
     source("scripts/m_GMC_bears.Rev")
-
+{% endsnippet %}
 We can now create our workspace model variable with our fully specified
 model DAG. We will do this with the `model()` function and provide a
 single node in the graph (`er`).
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
     mymodel = model(er)
-
+{% endsnippet %}
 *Run the Power-Posterior Sampler and Compute the Marginal Likelihoods* —
 With a fully specified model, we can set up the `powerPosterior()`
 analysis to create a file of 'powers' and likelihoods from which we can
@@ -414,29 +418,29 @@ sampling closer and closer to the prior as the power decreases.
 
 First, we initialize a monitor which will log the MCMC samples for each
 parameter at every step in the power posterior.
-
-    monitors[1] = mnModel(filename="output/GMC_posterior_pp.log",printgen=10, separator = TAB)
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
+    monitors.append( mnModel(filename="output/GMC_posterior_pp.log",printgen=10, separator = TAB) )
+{% endsnippet %}
 Next, we create the variable containing the power posterior. This
 requires us to provide a model and vector of moves, as well as an output
 file name. The `cats` argument sets the number of power steps. Once we
 have specified the options for our sampler, we can then start the run
 after a burn-in/tuning period.
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
     pow_p = powerPosterior(mymodel, moves, monitors, "output/GMC_bears_powp.out", cats=50, sampleFreq=10) 
     pow_p.burnin(generations=5000,tuningInterval=200)
     pow_p.run(generations=1000)  
-
+{% endsnippet %}
 Compute the marginal likelihood using two different methods,
 stepping-stone sampling and path sampling.
-
+{% snippet scripts/mlnl_GMC_bears.Rev %}
     ss = steppingStoneSampler(file="output/GMC_bears_powp.out", powerColumnName="power", likelihoodColumnName="likelihood")
     ss.marginal() 
 
     ### use path sampling to calculate marginal likelihoods
     ps = pathSampler(file="output/GMC_bears_powp.out", powerColumnName="power", likelihoodColumnName="likelihood")
     ps.marginal() 
-
+{% endsnippet %}
 If you have entered all of this directly in the RevBayes console, you
 will see the marginal likelihoods under each method printed to screen.
 Otherwise, if you have created the separate `Rev` file in the `scripts`

tutorials/clocks/scripts/m_GMC_bears.Rev

@@ -9,7 +9,7 @@
 ### the clock rate applied to every branch in the tree
 ### with a gamma prior and a scale move
 clock_rate ~ dnGamma(2.0,4.0)
-moves.append(mvScale(clock_rate,lambda=0.5,tune=true,weight=5.0))
+moves.append( mvScale(clock_rate,lambda=0.5,tune=true,weight=5.0) )
 
 ### set up the GTR model and instantaneous rate matrix from file
 source("scripts/m_GTR.Rev")

tutorials/clocks/scripts/m_GTR.Rev

@@ -8,17 +8,15 @@
 
 ### the stationary base frequencies 
 ### with a flat Dirichlet prior
-sf_hp <- v(1,1,1,1)
-sf ~ dnDirichlet(sf_hp)
+sf ~ dnDirichlet(v(1,1,1,1))
 
 ### the exchangeability rates 
 ### with a flat Dirichlet prior
-er_hp <- v(1,1,1,1,1,1)
-er ~ dnDirichlet(er_hp)
+er ~ dnDirichlet(v(1,1,1,1,1,1))
 
 ### a deterministic node for the instantaneous rate matrix
 Q := fnGTR(er,sf)
 
 ### simplex moves for er and sf
-moves.append(mvSimplexElementScale(er, alpha=10.0, tune=true, weight=3.0))
-moves.append(mvSimplexElementScale(sf, alpha=10.0, tune=true, weight=3.0))
+moves.append( mvSimplexElementScale(er, alpha=10.0, tune=true, weight=3.0) )
+moves.append( mvSimplexElementScale(sf, alpha=10.0, tune=true, weight=3.0) )

tutorials/clocks/scripts/mcmc_TimeTree_bears.Rev

@@ -31,9 +31,9 @@ source("scripts/m_UCLN_bears.Rev")
 mymodel = model(er)
 
 ### set up the monitors that will output parameter values to file and screen 
-monitors[1] = mnModel(filename="output/TimeTree_bears_mcmc.log", printgen=10)
-monitors[2] = mnFile(filename="output/TimeTree_bears_mcmc.trees", printgen=10, timetree)
-monitors[3] = mnScreen(printgen=100, root_time, tmrca_Ursidae)
+monitors.append( mnModel(filename="output/TimeTree_bears_mcmc.log", printgen=10) )
+monitors.append( mnFile(filename="output/TimeTree_bears_mcmc.trees", printgen=10, timetree) )
+monitors.append( mnScreen(printgen=100, root_time, tmrca_Ursidae) )
 
 ### workspace mcmc ###
 mymcmc = mcmc(mymodel, monitors, moves)

tutorials/clocks/scripts/mlnl_GMC_bears.Rev

@@ -22,7 +22,7 @@ source("scripts/m_GMC_bears.Rev")
 mymodel = model(er)
 
 ### create a monitor to log the samples to file
-monitors[1] = mnModel(filename="output/GMC_posterior_pp.log",printgen=10, separator = TAB)
+monitors.append( mnModel(filename="output/GMC_posterior_pp.log",printgen=10, separator = TAB) )
 
 ### compute power posterior distributions
 pow_p = powerPosterior(mymodel, moves, monitors, "output/GMC_bears_powp.out", cats=50, sampleFreq=10) 

tutorials/clocks/scripts/mlnl_UCLN_bears.Rev

@@ -22,7 +22,7 @@ source("scripts/m_UCLN_bears.Rev")
 mymodel = model(er)
 
 ### create a monitor to log the samples to file
-monitors[1] = mnModel(filename="output/UCLN_posterior_pp.log",printgen=10, separator = TAB)
+monitors.append( mnModel(filename="output/UCLN_posterior_pp.log",printgen=10, separator = TAB) )
 
 ### compute power posterior distributions
 pow_p = powerPosterior(mymodel, moves, monitors, "output/UCLN_bears_powp.out", cats=50, sampleFreq=10) 

tutorials/clocks/tests.txt

@@ -0,0 +1 @@
+mcmc_TimeTree_bears.Rev

tutorials/sse/scripts/mcmc_fBiSSE.Rev

@@ -32,6 +32,7 @@ root_age <- tree.rootAge()
 # vectors for moves and monitors
 moves    = VectorMoves()
 monitors = VectorMonitors()
+seed(3)
 
 ###
 # create parameters for the diversification and serial sampling rates