#
#
# 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^3

w_rbeta = function(nSim,alpha,beta){ 
  X = NULL
  while( length(X)<nSim ){
    nGen = nSim - length(X)
    U = runif(nGen)
    V = runif(nGen)
    UtoAI = U^(1/alpha)
    VtoBI = V^(1/beta)
    rat   = UtoAI / ( UtoAI + VtoBI ) 
    X     = c(X,rat[ UtoAI + VtoBI <= 1 ])
  }
  X = X[1:nSim]

  X
}

alpha = 2.1
beta = 3.1 

X = w_rbeta(nSim,alpha,beta)

###postscript("../../WriteUp/Graphics/Chapter2/ex_2_16_small_beta.eps", onefile=FALSE, horizontal=FALSE)

xRng = seq( min(X), max(X), length.out=100 )
hist(X,freq=F,breaks=100,main=sprintf("Small alpha= %3.1f, beta= %3.1f",alpha,beta))
lines(xRng,dbeta(xRng,alpha,beta),lwd=2,col="sienna")

###dev.off()


alpha = 7.1
beta = 9.1 

X = w_rbeta(nSim,alpha,beta)

###postscript("../../WriteUp/Graphics/Chapter2/ex_2_16_large_beta.eps", onefile=FALSE, horizontal=FALSE)

xRng = seq( min(X), max(X), length.out=100 )
hist(X,freq=F,breaks=100,main=sprintf("Large alpha= %3.1f, beta= %3.1f",alpha,beta))
lines(xRng,dbeta(xRng,alpha,beta),lwd=2,col="sienna")

###dev.off()