function(alpha, beta, n, s)
{
	lbeta <- function(a1, a2)
	{
		lgamma(a1) + lgamma(a2) - lgamma(a1 + a2)
	}
	lcon <- lgamma(n + 1) - lgamma(s + 1) - lgamma(n - s + 1)
	pred <- exp(lcon + lbeta(s + alpha, n - s + beta) - lbeta(alpha, beta))
	print(cbind(s, pred))
	X11()    # motif()
	plot(s, pred, type = "p", xlab = paste("number of successes in", n,
		"trials"), ylab = "predictive density")
}