# epage 74
#
# 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;

# Part (a): The code for the Poisson RV generator: 
#
lambda = 100
spread = 4*sqrt(lambda)
t = round( seq( max(0,lambda-spread),lambda+spread,1 ) )
prob = ppois(t,lambda)
X = rep(0,Nsim)
for( i in 1:Nsim ){
  u = runif(1)
  X[i] = t[1] + sum(prob<u)
}

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

xRng = seq(min(X),max(X),1)
brks = c( xRng[1] - 0.5, xRng + 0.5 )
hist(X,freq=F,breaks=brks,main="Poisson")

lines( xRng, dpois(xRng,lambda), col="blue")

###dev.off()



# Part (b): The negative binomial distribution:
#
mc_negBinomial = function(Nsim,r,p){
  # a function to generate negative binomial random variables using the inverse transform method
  # 

  nb_mean = r*(1-p)/p
  nb_var  = nb_mean/p

  spread = 4*sqrt(nb_var)
  t = round( seq( max(0,nb_mean-spread),nb_mean+spread,1 ) )

  w_negbin = choose( r + t - 1, t ) * p^r * (1-p)^t
  prob = cumsum( w_negbin ) 
  
  X = rep(0,Nsim)
  for( i in 1:Nsim ){
    u = runif(1)
    X[i] = t[1] + sum(prob<u)
  }
  
  X
}


r = 10


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

par(mfrow=c(1,3))

for( p in c(0.01,0.1,0.5) ){

  X = mc_negBinomial(Nsim,r,p)

  xRng = seq(min(X),max(X),1)
  brks = c( xRng[1] - 0.5, xRng + 0.5 )
  hist(X,freq=F,breaks=brks,main=sprintf("NB:r=%1d; p=%5.2f",r,p))

  w_negbin = choose( r + xRng - 1, xRng ) * p^r * (1-p)^xRng
  lines( xRng, w_negbin, col="blue" )
}
  
###dev.off()


# Part (c): The logarithmic series 
#
mc_logarithmic = function(Nsim,p){
  # a function to generate logarithmic random variables using the inverse transform method
  # 

  log_mean = ((1-p)/p)*(-1/log(p))
  log_var  = (-(1-p))*(1-p+log(p))/( p^2 * log(p)^2 )

  spread = 4*sqrt(log_var)
  t = round( seq( max(1,log_mean-spread),log_mean+spread,1 ) ) # the smallest value for t is now *1*

  w_log = ( - (1-p)^t )/( t * log(p) )
  prob = cumsum( w_log ) 
  
  X = rep(0,Nsim)
  for( i in 1:Nsim ){
    u = runif(1)
    X[i] = t[1] + sum(prob<u)
  }
  
  X
}

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

par(mfrow=c(1,3))

for( p in c(0.001,0.01,0.5) ){

  X = mc_logarithmic(Nsim,p)
  xRng = seq(min(X),max(X),1)
  brks = c( xRng[1] - 0.5, xRng + 0.5 )
  hist(X,freq=F,breaks=brks,main=sprintf("Log: p=%5.3f",p))
  
  w_log = ( - (1-p)^xRng )/( xRng * log(p) )
  lines( xRng, w_log, col="blue" )

}

###dev.off()