Example R programs and commands

	 30. Multinomial densities: calculation and sampling
    

# All lines preceded by the "#" character are my comments.
# All other left-justified lines are my input.
# All other indented lines are the R program output.

# MULTINOMIALS
#
# Setup:
#   r possible outcomes with probabilities p1,...,pr
#        -- so p1+...+pr = 1

#   n experiments yielding x1,...,xr outcomes
#        -- so x1+...+xr = n


# Evaluate the multinomial probability density function (pdf)
#
# Function:
#   dmultinom(c(x1,...,xr), prob=c(p1,...,pr))
#

# Compute single pdf values:
p1<- 0.2; p2<- 0.3; p3<- 0.5;    # 3 possible outcomes
dmultinom(c(5,5,5), prob=c(p1,p2,p3))  # prob. of 5 each
dmultinom(c(0,0,9), prob=c(p1,p2,p3))  # prob. of 9 3's
dmultinom(c(0,2,7), prob=c(p1,p2,p3))  # prob. of 2 2's and 7 3's.

# Sum with a loop to compute probability of no 1's in 9 trials:
p1<- 0.2; p2<- 0.3; p3<- 0.5; # 3 possible outcomes
ps <- c(p1,p2,p3)             # probabilities of outcomes 1,2,3
total <- 0;                   # initialize what gets the sum.  
for (k in 0:9)     # Loop and sum over (0,0,9),(0,1,8),...,(0,9,0)
   total <- total + dmultinom(c(0,k,9-k), prob=ps);
total



# Sample the multinomial probability density function (pdf)
#
# Function:
#   rmultinom( m, size=n,  prob=c(p1,...,pr))
#


# Generate sample runs of a 9-trial experiment:
p1<- 0.2; p2<- 0.3; p3<- 0.5;    # 3 possible outcomes
rmultinom(1, size=9, prob=c(p1,p2,p3))  # 1 run of 9 trials
rmultinom(5, size=9, prob=c(p1,p2,p3))  # 5 runs of 9 trials

# Estimate p1,p2,p3 by averaging over 50 runs:
sim <- rmultinom(50, size=9, prob=c(p1,p2,p3));  # 3x50 matrix
rowMeans(sim)       # average number of 1s, 2s, and 3s
rowMeans(sim)/9     # estimate for p1,p2,p3