#
#
# 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.
# 
#-----

# Part (a)
#

# the gamma:
# 
a = 2
beta = 2.4

Nsim = 10^5

T = matrix(rexp(a*Nsim,1),ncol=a)
X = beta * apply(T,1,sum) # these are our samples 

if( FALSE ){ 
  hist(X,freq=F,breaks=100)
  xRng = seq( min(X), max(X), length.out=100 )
  lines(xRng,dgamma(xRng,a,1/beta),lwd=2,col="sienna")
}

# the beta: 
#
b = 4

if( F ){ 
  # use different samples from Exp(1) for numerator and denominator ... this does NOT seem to be correct 
  #
  T1 = matrix(rexp(a*Nsim,1),ncol=a)
  num = apply(T1,1,sum)
  T2 = matrix(rexp((a+b)*Nsim,1),ncol=a+b)
  den = apply(T2,1,sum)
  X = num/den # these are our samples
}else{
  # use the *same* samples from an Exp(1) for the numerator and the denominator ... this seems to be correct 
  #
  T =  matrix(rexp((a+b)*Nsim,1),ncol=a+b) 
  num = apply(T[,1:a],1,sum)
  den = apply(T,1,sum)
  X = num/den # these are our samples
}
  

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

xRng = seq( min(X), max(X), length.out=100 )
hist(X,freq=F,breaks=100)
lines(xRng,dbeta(xRng,a,b),lwd=2,col="sienna")

###dev.off()