Orlando Schwery, 24. August 2021
This is based on the following tutorials:
If the RevBayes executable is on your system path, you should be able to open your command-line and start the application like this:
rb # or maybe ./rb depending on your systemNow we’re in the RevBayes terminal, just like we would be in R on command-line (or the vanilla R GUI), and we can type some simple calculations:
# Simple mathematical operators:
1 + 1 # Addition
10 - 5 # Subtraction
5 * 5 # Multiplication
10 / 2 # Division
2^3 # Exponentiation
5%2 # ModuloWe can also call some simple functions:
# Math functions
exp(1) # exponential function
ln(1) # logarithmic function with natural base
sqrt(16) # square root function
power(2,2) # power function: power(a,b) = a^bOften, for larger analyses, we probably won’t want to run RevBayes interactively, but save all our code in a script and execute that whenever we need to.
We can use the familiar functions from R to see where our working directory currently is, and can set a new one:
getwd()
setwd("/mnt/c/Users/oschwery/Documents/Teaching/RevBayes_Workshop_UI2021/")(If you’re on Windows, depending on the details, this might not work for you, for… reasons… We can look at this, but might have to go with some workarounds for now.)
Once you’re in the right directory, you can run a script therein by sourcing it:
source("testcode1.Rev")Once you’re done, you can leave RevBayes like this:
q()You can also call the code file directly from the terminal as you call the RevBayes program:
rb /mnt/c/Users/oschwery/Documents/Teaching/RevBayes_Workshop_UI2021/testcode2.RevNote that for this one, RevBayes quit by itself, because we added the q() as the last line of the code. Also note that we had to provide the whole path to that file here, because if we call RevBayes like that, the working directory will still be at the default location and not where we’ve set it in the previous session (some OIs might automatically set it to the directory from which rb was called though).
First off, let’s focus on the different ways to assign variables, since the types of variables are crucial here.
In R, we are used to two ways to assign, and they are largely interchangeable:
x <- 1
y = 2In Rev, we have a different assigner for each of the different nodes:
p <- 0.5 # constant variable
x ~ dnBernoulli(p) # stochastic variable
x := exp(p) # deterministic variable
x.clamp(1) # clamped variable
y = 2 # inference (i.e. non-model) variable
for (i in 1:N){ x[i] ~ dnBernoulli(p)} # plateIf we tried to run x.clamp(1) above, we would be getting an error stating that it doesn’t know that function. The reason for this is that we’ve overwritten the stochastic variable x with the deterministic variable x. We’ll discuss the logic of clamping below, but the TLDR is that we assign observed data to an object, which doesn’t really make sense if it’s a deterministic variable (as the name suggests, it is deterministic, so exp(p) will always be the same as long as p stays the same, so there’s no sense in attaching an observed value to this).
Rev telling us that it doesn’t know the function clamp is unfortunate for us trying to understand the problem, but reflects the fact that functions are tied to the types of nodes they can be applied to ( i.e. Rev is an object-oriented language), and Rev is essentially telling us “I don’t know a function called clamp that I could apply to this x here (which is a deterministic node/variable)”.
If we re-run x as the stochastic variable and try to clamp it now, it works:
x ~ dnBernoulli(p) # stochastic variable
x.clamp(1) # clamped variableNote that as in R, functions are written in the format funcname(args), but unlike R, if a function is applied to an object, it isn’t funcname(obj, args), but obj.funcname(args).
As one would expect, if we assign a stochastic variable repeatedly, the value changes every time:
for( i in 1:12 ){
e[i] ~ dnExponential(1.0)
}
eBut also note how unlike in R, we didn’t have to initiate the variable e before starting the for-loop.
For loops are a very useful programming tool to do the same operation for a large number of objects:
# abstract structure
for ([variable] in [set of values]) {
[statements using variable]
}
# an example that will simply print the current variable i (for iteration) in each round:
for (i in 1:10) {
i
}
# an example that will add a new value to each of two vectors in each round:
lambda <- 2
for (i in 1:10) {
y[i] ~ dnUniform(0, 1)
z[i] := -ln(y[i]) / lambda
}This last example will throw an error, because
1) there already IS an object y that we created above, and
2) that object was assigned with a =, and not a ~, which is why Rev can’t overwrite it.
This is a fundamental difference between R and Rev, because based on how a variable was assigned, Rev will keep track of what type of node/variable that was, and will prevent it from being treated like a different type of node. So here it won’t let us assign draws from the uniform distribution to y, because it isn’t a stochastic node, but an inference variable (or as it calls it a variable of type 'Natural'). In R, the old y would have just been overwritten anyways.
Also, to clarify, if we had just tried to overwrite this old y with a new one like y ~ dnUniform(0, 1), that would have worked. But instead we tried to add to an index of the existing one (y[i]), and that’s what we couldn’t do because of the different variable types.
Now let’s try this again, renaming this y to Y:
lambda <- 2
for (i in 1:10) {
Y[i] ~ dnUniform(0, 1)
z[i] := -ln(Y[i]) / lambda
}Again, we didn’t need to initalise the objects Y and z.
Above we renamed the new y to Y to avoid conflict, but we could also have deleted it using clear(). You can empty the environment with clear() and delete a particular object with clear(a), although it can have some weird behaviour sometimes:
ls() # shows us the variables we just had
clear() # removes them all
ls() # as we can see...
# set a few new objects
a <- 1
b = 2
c := a+b
ls() # see them
# clear each individually
clear(a)
clear(b)
clear(c)
# depending on your operating system, you may have gotten an error message when trying to clear b...
ls() # ...which is now still around...
clear() # clearing everything should take care of that though
ls() # are they gone?We can now attempt to code up the lognormal model we’ve seen on the PowerPoint (slide 10), starting with the gamma distribution on the top left, giving us M:
alpha <- 2.0
beta <- 4.0
M ~ dnGamma(alpha, beta)Followed by the exponential distribution on the right, giving us sigma:
lambda <- 1.0
sigma ~ dnExponential(lambda)For both of these, we created constant nodes (<-) which we then used as parameters for the stochastic nodes (~) that let us draw from these distributions. Now finally, we use both M and sigma to create the deterministic node mu, which is just a simple calculation of the former two:
mu := ln(M) - (sigma^2)/2.0Now we could use mu and sigma to draw xs from a lognormal distribution like this: x ~ dnLognormal(mu, sigma). However, we don’t want to use this model to draw variables, but to evaluate data.
Let’s use this as an attempt to exemplify clamped nodes. First we create a vector of seven values drawn from that lognormal distribution (which is conincidentally the same amount as our data).
for( i in 1:7 ) {
x[i] ~ dnLognormal(mu, sigma)
}
x # inspecting xIf we look at x now, we see seven random values drawn from that distribution. But now we got seven observations here instead that we’re clamping to this:
observations <- [0.20, 0.21, 0.03, 0.40, 0.65, 0.87, 0.22] # note this is one way to create a vector
N <- observations.size() # size() is a similar function to R's length()
# it is good practice to code the loop such that
# it automatically goes through as many iterations
# as the target object has elements
# (e.g. vector entries here), so we don't have to
# keep changing the code in case our data changes
for( i in 1:N ){
x[i].clamp(observations[i])
}
x # inspecting how x has changedLooking at x again, we see that it now just has the values of our observations. The logic here is the following:
The model we created relies on the underlying distribution, which determines with what probability we draw a particular number from it. Now we might want to turn this around and evaluate for a given set of numbers what the probability is that they would be drawn from that distribution. By clamping those values to the object with the drawn numbers, we’re essentially connecting them with the distribution, i.e. we’re telling the model what values were drawn, instead asking it to draw them, if you will. [If you let that sink in, it makes sense that we can probably only clamp observed data to stochastic nodes…]
This connects to the points above, where we said that Rev keeps tabs on what type a variable is, and only allows them to be manupulated in accordance with that type (see clamping to a deterministic variable xor adding stocastic values to a non-model y above). So when we look at the observations and x now, they may look the same:
observations
xBut they are not, as one is a vector of constant variables and the other is clamped/stochastic. This can be shown by trying to get other information from it, e.g. the probability of a particular value under the that lognormal distrubution we based this on:
observations[1]
x[1]
observations[1].probability()
x[1].probability()Again, Rev claims to not know probability when applied to the constant node, but gives us a value when applied to the clamped one.
We can also apply probability to stochastic nodes without observations clamped to them:
lambda <- 1.0 # assign constant node `lambda' with value `1'
x ~ dnExponential(lambda) # create stochastic node with exponential distribution and parameter `lambda'
x # print value of stochastic node `x'
x.probability() # print the probability if `x'
x.lnProbability() # print the log-probability if `x'Finally, a point about the modularity of RevBayes: We could now exchange the M above with a binomial lognormal distribution instead:
mean_1 <- 0.5
mean_2 <- 2.0
sd_1 <- 1.0
sd_2 <- 1.0
weight <- 0.5
M ~ dnBimodalLognormal(mean_1, mean_2, sd_1, sd_2, weight)Now we could just rerun the above lines giving us mu and the xs again unchanged and would get new results based on a different mean. This is trivial for a small model like this, but can become really handy later on.
Just like R, Rev is case-sensitive:
exp(1) # correct lower case name
Exp(1) # wrong upper case nameVectors are created and manipulated as follows:
z <- v(1.0, 2.0, 3.0) # create a vector (but see how we created it in another way above)
z[1] # print the first entry
z[1] <- 10 # change the value of the first entry
zMultiple statements can go on one line if separated with a semicolon:
1 + 1; 2 + 2 # Multiple statements in one lineAnd as you probably noticed, coments can be entered with # and Rev will ignore them, but careful:
# try this in R
rawr <- c( 1,2,3, #some comment here
4,5,6,7)# try this in Rev
rawr <- [1,2,3, #some comment here
4,5,6,7]Another thing to be aware of is that there is no lazy-evaluation in Rev, unlike in R. This means that even if a section of the code shouldn’t be evaluated (e.g. because the condition for it to run is currently FALSE), Rev will still read it and potentially throw an error if that part of the code isnt’in order (which can really keep you searching for the problem if you’re unaware). [I meant to have an example here but couldn’t recreate such a case, so you’ll have to take my word for it…]
Just as a last tidbit of I ran into that can be hard to figure out based on the error message (or lack thereof): Apparently, there is a maximum length for how much can go into the same function. So if you want to e.g. write a whole range of calculations directly to file, your code may crash for mysterious reasons…
You can trigger help topics by adding a question mark, just as in R:
?dnNorm
?mcmc
?mcmc.runHowever, there’s limitations on how easily help can be gotten when compared to R…
Additionally, if you type a function without brackets, Rev will show you the proper usage of it:
mcmcAlso see RevBayes User’s Forum for more help.