# The function "prob.mass(x,prob)" draws the actual probability mass
# function on the top of histogram.
prob.mass = function(x,prob,...){
  for(i in 1:length(x)){
    polygon(c(x[i]-0.5,x[i]-0.5,x[i]+0.5,x[i]+0.5),
            c(0,prob[i],prob[i],0),...)
  }
}
# If the distribution is binomial, a vector of probabilities is
# automatically generated as in the following example. 
n = 20
p = 0.1
x = 0:n
prob = dbinom(x,n,p)
barplot(prob,names.arg=x,col=2)
mu = sum(x * prob)
cat("\n E[X] =", mu)
sigma.square = sum((x - mu)^2 * prob)
cat("\n Var(X) =", sigma.square)
# Also, it becomes much simpler to generate a random sample for the
# binomial distribution.
size = 200
sample = rbinom(size,n,p)
cat("\n sample =", sample)
hist(sample, freq=F, breaks=seq(-0.5,n+0.5,by=1.0), ylim=c(0,max(prob)+0.1), col='green')
prob.mass(x,prob,border='red',lwd=2)
cat("\n sample mean =", mean(sample))
cat("\n sample var =", var(sample))