ST 321 Simulate Discrete Random Variable March 26, 2004

In this lab session you will use the Inverse Transform method to simulate different discrete random variables.

Task 1: In one of our previous classes we have studied the following R/Splus function:

inverse.transform.discrete=function(x,p) {

#

# generate discrete random variables using the

# inverse transformation method

#

# x: possible values for the random variable

# p: probabilities

#

# x and p must be vectors of the same length

#

psofar=p[1]

j=1

X=x[1]

U=runif(1)

while (U>psofar) {

j=j+1

psofar=psofar+p[j]

X=x[j]

}

return(X)

}

Use the above R function to simulate 10 random variables with the probability density function: P(X=1)=0.2, P(X=2)=0.1, P(X=3)=0.35, P(X=4)=0.2 and P(X=5)=0.15. (Note: you may want to consult the corresponding handout from class.)

Task 2: In class we have shown that geometric random variables can be simulated by one single R/Splus statement. Verify this by running the following R/Splus commands.

(Note: the following commands first define a function my.rgeom for generating geometric random variables. Then it simulates n=1000 geometric random variables (with success probability 0.3) from my.rgeom and another n=1000 from the R built-in function rgeom. Finally it compares the two sets of simulated values by plotting their histograms.)

my.rgeom=function(n,prob) {

return(floor(log(runif(n))/log(1-prob))+1)

}

n=1000

p=0.3

r.values=rgeom(n,p)+1

# note: in R rgeom uses a different parametrization from ours

# we need to add 1 to make them the same

my.values=my.rgeom(n,p)

par(mfrow=c(2,1)) # to split the graphic window into two

hist(r.values,probability=T,nclass=20,ylim=c(0,0.6),xlim=c(1,21))

hist(my.values,probability=T,nclass=20,ylim=c(0,0.6),xlim=c(1,21))

Task 3: In class we have also discussed the generation of Poisson random variables. The key observation behind our generation method is the identity P(X=i+1) = P(X=i)*lambda/(i+1) . The following R/Splus function implement our method for generating Poisson random variables:

my.rpois=function(lambda) {

# generate 1 poisson random variable

U=runif(1)

i=0

probi=exp(-lambda)

probsum=probi

while (U > probsum) {

probi=probi*lambda/(i+1)

probsum=probsum+probi

i=i+1

}

return (i)

}

Now run the following commands:

n=1000

lambda=3.5

r.values=rpois(n,lambda)

my.values=seq(0,0,length=n)

for (i in 1:n) {

my.values[i]=my.rpois(lambda)

}

hist(r.values,probability=T,nclass=10,ylim=c(0,0.3),xlim=c(0,10))

hist(my.values,probability=T,nclass=10,ylim=c(0,0.3),xlim=c(0,10))

Task 4: Modify my.rpois to generate binomial random variables. Use your modified function to generate 1000 binomial random variables with N=10 trials and success probability p=0.7. Compare your simulated values with those from the built-in R function rbinom as before (i.e., by plotting histograms).

Bonus Task: How would you generate hypergeometric and negative binomial random variables? Hint: find the relation between P(X=i) and P(X=i+1).