#
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

Nsim = 10^4; n = 25; p = 0.2;


wrbinom = function(Nsim,n,p){
  # a function to generate binomial random variables using the inverse transform method
  # 
  
  prob = pbinom(0:n,n,p)
  X = rep(0,Nsim)
  for( i in 1:Nsim ){
    u = runif(1)
    X[i] = 0 + sum(prob<u)
    }

  X
}

system.time(wrbinom(Nsim,n,p)) # My generator 

system.time(rbinom(Nsim,n,p)) # The R generator is rbinom


# Part (b):
#

wUniformAlpha = function(Nsim,alpha){
  # a function to generate random variables from U_{[0,alpha]}
  # 

  X = rep(0,Nsim)
  for( i in 1:Nsim ){
    u = runif(1)
    while( u > alpha ) u = runif(1)
    X[i] = u 
  }

  X
}

# plot the histogram of these data points  
X = wUniformAlpha(Nsim,0.2)
hist(X,freq=F)

# compare the time to simulate 
Nsim = 10^5

alpha = 0.1
alpha = 0.9

system.time( wUniformAlpha(Nsim,alpha) )

system.time( runif(Nsim,max=alpha) )

system.time( alpha*runif(Nsim) )