# 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),...)
  }
}
# The expectation calculation in theory.
# Given the freqency function
# p(1) = 0.2; p(2) = 0.3; p(3) = 0.4; p(4) = 0.1
# we can calculate E[X], Var(X), and SD of X.
prob = c(0.2, 0.3, 0.4, 0.1)
x = c(1,2,3,4)
barplot(prob, names.arg=x, col=2)
x * prob
mu = sum(x * prob)
cat("\n E[X] =", mu)
sigma.square = sum((x - mu)^2 * prob)
cat("\n Var(X) =", sigma.square)
cat("\n SD =", sqrt(sigma.square))
# Now we generate a sample of size 200,
# and calculate "average", "sample variance", and "standard deviation"
size = 200
# sample <- as.numeric(x %*% rmultinom(size,1,prob))
sample = sample(x, size=size, replace=T, prob=prob)
cat("\n sample =", sample)
hist(sample, freq=F, breaks=seq(min(x)-0.5,max(x)+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))
cat("\n sample SD =", sd(sample))