### plot Bin(12, 1/3), like Weldonâ€™s dice

x <- seq(0,12)
plot(x=x, y=dbinom(x,12,1/3), type="h", lwd=8)

better labels, and lines()

x <- seq(0,12)
plot(x=x, y=dbinom(x,12,1/3), type="h",
xlab="Number of successes", ylab="Probability",
lwd=8, lend="square")
lines(x=x, y=dbinom(x,12,1/3))

use Weldonâ€™s dice data, calculate expected frequencies

data(WeldonDice, package="vcd")
Weldon.df <- as.data.frame(WeldonDice)   # convert to data frame

Prob <- dbinom(x, 12, 1/3)               # binomial probabilities
Prob <- c(Prob[1:10], sum(Prob[11:13]))  # sum values for 10+
Exp= round(sum(WeldonDice)*Prob)         # expected frequencies
Diff = Weldon.df[,"Freq"] - Exp          # raw residuals
Chisq = Diff^2 /Exp
data.frame(Weldon.df, Prob=round(Prob,6), Exp, Diff, Chisq)
##    n56 Freq     Prob  Exp Diff     Chisq
## 1    0  185 0.007707  203  -18 1.5960591
## 2    1 1149 0.046244 1216  -67 3.6916118
## 3    2 3265 0.127171 3345  -80 1.9133034
## 4    3 5475 0.211952 5576 -101 1.8294476
## 5    4 6114 0.238446 6273 -159 4.0301291
## 6    5 5194 0.190757 5018  176 6.1729773
## 7    6 3067 0.111275 2927  140 6.6962761
## 8    7 1331 0.047689 1255   76 4.6023904
## 9    8  403 0.014903  392   11 0.3086735
## 10   9  105 0.003312   87   18 3.7241379
## 11  10   18 0.000544   14    4 1.1428571

### create a set of binomial distributions, Bin(12, p)

# dbinom can take vector inputs for the parameters

x <- seq(0,12)                # values of x
p <- c(1/6, 1/3, 1/2, 2/3)    # values of p
x <- rep(x, 4)                # replicate, length(p) times
p <- rep(p, each=13)          # replicate, each length(x)
bin.df <- data.frame(x, prob = dbinom(x, 12, p), p)
bin.df$p <- factor(bin.df$p, labels=c("1/6", "1/3", "1/2", "2/3"))
str(bin.df)
## 'data.frame':    52 obs. of  3 variables:
##  $x : int 0 1 2 3 4 5 6 7 8 9 ... ##$ prob: num  0.1122 0.2692 0.2961 0.1974 0.0888 ...
##  $p : Factor w/ 4 levels "1/6","1/3","1/2",..: 1 1 1 1 1 1 1 1 1 1 ... best version: using expand.grid() XP <-expand.grid(x=0:12, p=c(1/6, 1/3, 1/2, 2/3)) bin.df <- data.frame(XP, prob=dbinom(XP[,"x"], 12, XP[,"p"])) bin.df$p <- factor(bin.df$p, labels=c("1/6", "1/3", "1/2", "2/3")) str(bin.df) ## 'data.frame': 52 obs. of 3 variables: ##$ x   : int  0 1 2 3 4 5 6 7 8 9 ...
##  $p : Factor w/ 4 levels "1/6","1/3","1/2",..: 1 1 1 1 1 1 1 1 1 1 ... ##$ prob: num  0.1122 0.2692 0.2961 0.1974 0.0888 ...

### Using lattice::xyplot() to plot the distributions

library(lattice)

mycol <- palette()[2:5]
xyplot( prob ~ x, data=bin.df, groups=p,
xlab=list('Number of successes', cex=1.25),
ylab=list('Probability',  cex=1.25),
type='b', pch=15:17, lwd=2, cex=1.25, col=mycol,
key = list(
title = 'Pr(success)',
points = list(pch=15:17, col=mycol, cex=1.25),
lines = list(lwd=2, col=mycol),
text = list(levels(bin.df\$p)),
x=0.9, y=0.98, corner=c(x=1, y=1)
)
)