gammasim <- 
function(r, n)
{
#
# GAMMASIM simulates n independent realisations of a gamma(r,1) random
#   variable, using rejection.
#
	y <- rep(0, n)
	for(j in 1:n) {
		t <- 2	# inefficient SPLUS programming
		u <- 3	# values to start the 'while' command
		while(t < u) {
			e <-  - r * log(runif(1))
			t <- ((e/r)^(r - 1)) * exp((1 - r) * (e/r - 1))
			u <- runif(1)
		}
		y[j] <- e
	}
	y
}
binsim <-
function(m, p, n)
{
#
# BINSIM  simulates n independent realisations of a Bin(m,p) random
#   variable.
#
	U <- matrix(runif(m * n), nrow = m)
	Uindx <- U - p < 0
	apply(Uindx,2,sum)
}
betasim <- 
function(a, b, n)
{
#
# BETASIM simulates n independent realisations of a Be(a,b) random
#   variable.
#   Calls the function GAMMASIM.
#
	x1 <- gammasim(a, n)
	x2 <- gammasim(b, n)
	x1/(x1 + x2)
}
# Program to perform Gibbs sampling for Example 7.1, and fixed n.
# 'm': denotes the number of cycles before sampling
# 'k': denotes the number of repeats to be used
# 'n','alpha' and 'beta' are the known parameters of the
# bivariate distribution.
# Calls the functions BETASIM and BINSIM.
#_______________________________________________________________________
	m <- 50
	k <- 100
	n <- 20
	alpha <- 1
	beta <- 2
	y1 <- rep(0, k)
	x1 <- rep(0, k)
	for(j in 1:k) {
		x <- n * runif(1)	# start for x
		for(i in 2:m) {
			y <- betasim(alpha + x, n - x + beta, 1)
			x <- binsim(n, y, 1)
		}
		y1[j] <- y	# we record the values
		x1[j] <- x	# after m cycles
	}
	X11()
	hist(y1, xlab = "y")
	X11()
	hist(x1, xlab = "x")