ST 321 Simulating
Experiments
Computing in Weber Classrooms 205-206:
1.
To log in, make
sure that the DOMAIN NAME is set to MATHSTAT.
2.
Use the class
username: st321. The password was distributed in the last lab session.
3.
On most of the
machines, your work should be saved to the folder C:\temp (on some machines,
this folder will be D:\temp). This is a temporary location for class work. The folder is emptied periodically. Work that you would like to keep should be
stored on a diskette or a memory stick.
4.
Course material
(data, homework files, etc.) will be placed in a subfolder st321/Spring2004 on
drive G:
5.
Folders can be
examined and navigated using Explore (click the right button over the Start
button and select Explore).
6.
At the end of class,
make sure that you logout from the computer (Start > Shutdown > Close all
programs and logon as a different user.)
Splus/R:
In this lab session you will see how to do comparisons
in Splus/R and how to write if-else statements. You will also see how to use Splus/R to simulate real experiments. Finally you will see how to estimate
probabilities using simulated repeated experiments. First, try out the following Splus/R comparison commands:
Comparisons:
>
set.seed(321)
>
x <- sample(1:6,20,replace=T) [simulate “roll a die 20 times”]
>
sum(x==2) [count
number of “two”’s]
>
sum(x!=1) [count
number of “non-one”’s]
>
20-sum(x==1) [same as
above]
>
sum(x>4) [count
number of “five”’s and “six”’s]
>
sum(x>=5) [same
as above]
>
max(x) [maximum
of x]
>
min(x) [minimum
of x]
>
T&F [&:
logical AND]
>
T&T
>
T|F [|:
logical OR]
>
F|F
>
#
>
#this is a comment [everything
after “#” will be ignored]
>
#ignored by Splus/R
>
#
Random exponentials:
>
x <- rexp(10, 0.5) [generate
10 exponential random]
>
x [variables
with rate 0.5]
If-Else Statements:
The following Splus/R program simulates the following experiment. First roll a fair die. If face one comes up, toss a fair coin one
time, if faces two or three comes up, toss a fair coin two times, otherwise
toss a fair coin three times. You can
copy and paste this program from the subfolder st321/Spring2004 on G: drive,
filename “labsession2”. (To run the
program, type experiment1()in Splus/R.)
experiment1
<- function()
{
face <-
sample(1:6,1,replace=T)
if (face==1) {
coin <-
sample(0:1,1,replace=T) # 0 for Tail, 1
for Head
} else if ((face==2)|(face==3))
{
coin <-
sample(0:1,2,replace=T)
} else {
coin <-
sample(0:1,3,replace=T)
}
numberH
<- sum(coin==1)
return(numberH)
}
Experiment 1 (cont):
Do you know how to calculate
P(2H), the probability of getting 2 heads? The following program simulates repeated
experiments to estimate this probability. Again, you can find this program in
the file “labsession2”. Note the cat command is
for displaying results on the screen.
estimate1
<- function(N=1000)
{
zerohead <- 0
onehead <- 0
twoheads <- 0
threeheads
<- 0
for (i in 1:N) {
x <- experiment1()
if (x==0) {
zerohead
<- zerohead + 1
} else if (x==1) {
onehead <- onehead + 1
} else if (x==2) {
twoheads
<- twoheads + 1
} else {
threeheads
<- threeheads + 1
}
}
cat("prob(0 H):", zerohead/N,
"\n")
cat("prob(1 H):", onehead/N,
"\n")
cat("prob(2 H):", twoheads/N,
"\n")
cat("prob(3 H):", threeheads/N,
"\n")
}
Task 1. Run estimate1 three times with three different values of N: 1000, 10000 and 100000. Compare the simulated values to the exact values: P(0H)=11/48, P(1H)=7/16, P(2H)=13/48 and P(3H)=1/16.
Task 2. The following program can be found in the file “labsession2”. What does it do?
estimate2
<- function(N=1000, rate=1, cutpoint=1)
{
x <- rexp(N,rate)
smaller <- x<cutpoint
return(sum(smaller)/N)
}
Task 3. Write a Splus/R program to generate a random variable X defined by the following procedure. First a fair coin is tossed 4 times. Let N be the total number of heads obtained. If N<2 (i.e., 0 or 1), throw a fair die; otherwise, throw a pair of fair dice. Let Y be the number obtained if one die is threw, or be the maximum of the two numbers obtained if a pair of dice are threw. Note that Y is always an integer between 1 and 6. Finally generate Y exponential random variables with rate 2. The random variable X is defined as the sum of these Y exponential random variables.
Task 4. Use the above
program to estimate P(X>3).