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