STAT3405: Computer Laboratory 1
STAT4066: Computer Laboratory 1
In this lab we will familiarise ourselves with the program OpenBUGS by doing some simulation
studies. The material in this lab follows closely material in the initial two chapters of Lunn
et al. (2012).
We will also learn how to run these BUGS programs from R by using either the BRugs package
or the rjags package. We will also start to learn how to use the probabilistic programming
language Stan.
Getting Started:
Log in at one of the PCs and start up the software package OpenBUGS. With Windows 10 on
the lab PCs, it is easiest to start a software package by typing its name into the “Type here
to search” field next to the Start button (bottom left of screen). Start typing OpenBUGS or
RStudio into this field, and select the appropriate app for launching as soon as it is found, to
start OpenBUGS or RStudio, respectively.
After starting OpenBUGS, you should see a window such as the one in Figure 1 appearing on
your screen.
Figure 1: OpenBUGS opening screen
Exercise 1: Coins: a Monte Carlo approach to estimating tail areas
Here we discuss how to solve Example 1 (taken from Chapter 1 of Lunn et al., 2012) given the
background reading material using OpenBUGS, following closely the discussion in Chapter 2
of Lunn et al. (2012). The model in that example was
Y ∼ Bin(0.5, 8)
1
and we want to know P (Y ≤ 2). This model is represented in the BUGS language as:
model {
Y ∼ dbin (0 .5 , 8)
P2 ← s tep (2 . 5 −Y ) # does Y = 2 , 1 , or 0?
}
P2 is a step function that will take on the value 1 if 2.5− Y ≥ 0, i.e., if Y is 2 or less, and 0 if
Y is 3 or more; this corresponds to the indicator function used in Example 1.
The following steps are used to run a basic model interactively in OpenBUGS. This process
will be explained in more details in later lectures/computer labs.
1. Make a new document (New from the File menu) and type in the BUGS model code.
2. Open Specification Tool from the Model menu. A dialog box like the one in Figure 2
will appear.
Figure 2: OpenBUGS Model Specification Tool
3. Highlight the word model in the BUGS code by double clicking on the word. Click on
check model in the model specification tool. Any error messages are shown on the
bottom left of the screen, or model is syntactically correct will appear if there
are no errors.
4. There are no observed data in this model; therefore we can ignore load data.
5. In this example it is sufficient to leave the number of parallel chains to run (num of
chains) at 1.
6. Click on compile. Again check for any error messages at the bottom left.
7. We can ignore load inits in this simple example.
8. Click on gen inits. A message initial values generated, models initialized
should appear (in the bottom left of the screen).
9. Open Update... from the Model menu (left of Figure 3) and Samples... from the
Inference menu (right of Figure 3)
10. Specify the nodes we want to monitor or record the sampled values for. In this case,
type P2 into the node box in the Sample Monitor Tool, and click set. Similarly, type
Y and click set.
2
Figure 3: OpenBUGS Update (left) and Sample Monitor (right) Tool
11. Type * into the node box in the Sample Monitor Tool, which means “all monitored
nodes”, and click trace to open a window where the sampled values will appear as they
are generated. (Note, this will not work in the current version of OpenBUGS, which
requires at least one update [next step] to have been performed before opening the trace
window.)
12. Go to the Update Tool and type the number of samples to be generated in the updates.
10,000 are sufficient in this example. Thus
(a) Either accept the default value of generating 1,000 values in each update step and
press update 10 times, or
(b) Enter the value 10,000 in the updates field and press update once.
To see traces of the generated values, you should select the first approach.
13. Type * in the Sample Monitor Tool again. Click stats to see summary statistics for
all monitored nodes (top of Figure 4), and density to see plots of empirical distributions
(bottom of Figure 4).
Figure 4: Summary statistics and empirical distributions for Y and P2 based on 10,000 simu-
lations.
Taking the empirical mean of P2 gives the estimated probability that Y will be 2 or fewer.
The summary statistics provided by OpenBUGS are shown on the top of Figure 4. The
mean and sd are simply the empirical average and standard deviation of the sample values
while, as described in the background material, the MC error (Monte Carlo error) provides an
assessment of the sampling error on the mean attributable to the limited number of iterations
performed. we note that the MC error calculated for P2 matches that obtained in Example 1 in
3
the background material. The 2.5%, median, and 97.5% values are the empirical percentiles,
while start is the iteration at which monitoring began, and sample indicates the total number
of iterations contributing to the summary statistics.
Notes on the BUGS language
The example above illustrates a number of aspects of the BUGS syntax. First, the entire
model description is enclosed in model{...}. Second, there are two types of “connectives”:
• <- represents a logical dependence. The left-hand side of a logical statement comprises
a logical node, and the right-hand side comprises an expression formed from the built-
in logical functions of OpenBUGS (see, among other, Lunn et al., 2012, Appendix B)
applied to a set of stochastic or logical notes, e.g. m <- a + b*x.
• ∼ represents stochastic dependence. The left-hand side of a stochastic statement com-
prises a stochastic node, and the right-hand side comprises a built-in distribution of
OpenBUGS (see, among other, Lunn et al., 2012, Appendix C), e.g., r ∼ dunif(a,b)
for a variable r that is uniformly distributed between a and b.
Note that in OpenBUGS, logical expressions are not permitted as parameters of distri-
butions, so a statement such as r ∼ dunif(2*a, b) is not permitted. However, such
an expression is allowed in JAGS.
• # is a comment character used to annotate the modelling code. Everything after a #
on the same line is ignored by BUGS. Clear and concise comments can be helpful when
reading and maintaining models, particularly if it is not immediately clear what a piece
of code does.
In general, each node in a model (apart from constants) should appear once and only once on
the left-hand side of a statement (although there are exceptions to this rule).
Exercise 2: Repairs: the “how many” trick
Suppose costs of a repair have a gamma distribution with a mean of A$100 and standard
deviation A$50. How many items will I be able to repair for A$1000? The gamma distribution
is a distribution with two parameters, typically denoted by α and β (and we will learn more
about this distribution later), but it is possible to work out that a gamma distribution with
parameters α = 4 and β = 0.04 has mean 100 and standard deviation 50.
The “how many trick” is then, at each iteration, to simulate costs Yi, i = 1, . . . , I, from this
distribution for sufficiently large I, find the empirical cumulative distribution, and the find
the value M which is the highest i such that the total cost does not exceed A$1000. We do
this by creating a new vector cum.step = 1, 2, . . . ,M, 0, . . . , 0, where M is the largest integer
such that
∑M
i=1 Yi < 1000, and use the ranked function to find the maximum of the elements
of cum.step:
model {
for ( i in 1 : 2 0 ) {
Y[ i ] ∼ dgamma(4 , 0 . 04 )
4
}
cum [ 1 ] ← Y[ 1 ]
for ( i in 2 : 2 0 ) {
cum [ i ] ← cum [ i −1] + Y[ i ]
}
for ( i in 1 : 2 0 ) {
cum.step [ i ] ← i ∗ s tep (1000 − cum [ i ] )
}
number ← ranked ( cum.step [ ] , 20) # maximum number in cum.step
check ← equa l s ( cum.step [ 2 0 ] , 0) # always 1 i f I=20 b i g enough
}
Running 10,000 iterations in OpenBUGS produces the summary statistics shown at the top of
Figure 5; the empirical distributions are shown at the bottom of that figure.
Figure 5: Summary statistics and empirical distributions of number of items which can be
repaired for A$1000, given a random item report cost of Y[1]. Also displays check that
simulation had large enough value of I.
Therefore the median number we can repair is 10, with a 95% predictive interval 7 to 13.
Note that Y (as demonstrated by monitoring Y[1]) has mean 100 and standard deviation 50,
as required. And I = 20 was large enough as the monitoring of check shows.
Exercise 3: BRugs
Download the R notebook Lab01-BRugs.Rmd from LMS (Content→ Worksheets for Computer
Practicals → Computer Lab 01) and open it from RStudio. Work through the notebook.
You will see that running a BUGS program using R and the BRugs package is very similar to
running the code in OpenBUGS. However, instead of issuing commands for the various steps
via the graphical user interface (GUI), we now issue corresponding commands in R.
Exercise 4: rjags
Download the R notebook Lab01-rjags.Rmd from LMS (Content→ Worksheets for Computer
Practicals → Computer Lab 01) and open it from RStudio. Work through the notebook.
5
Exercise 5: Stan
Download the R notebook Lab01-Stan.Rmd from LMS (Content→ Worksheets for Computer
Practicals → Computer Lab 01) and open it from RStudio. Work through the notebook.
Finishing Off:
When you’ve finished, make sure you first save all the work you want to keep and then close
down OpenBUGS by selecting Exit from the File menu. Remember to log off from your
computer before leaving.
References
Lunn, D., Jackson, C., Best, N., Thomas, A. and Spiegelhalter, D. (2012). The BUGS Book:
A Practical Introduction to Bayesian Analysis, Texts in Statistical Science, Chapman &
Hall/CRC, Boca Raton.
6































































































































































学霸联盟