xuebaunion@vip.163.com

3551 Trousdale Rkwy, University Park, Los Angeles, CA

留学生论文指导和课程辅导

无忧GPA：https://www.essaygpa.com

工作时间：全年无休-早上8点到凌晨3点

微信客服：xiaoxionga100

微信客服：ITCS521

openbugs/r代写-STAT3405

时间：2021-03-14

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

学霸联盟

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

学霸联盟