# set some initial parameters for the graphs
thebreaks <- seq(-1.5, 20.5, 1)
xlimits <- c(-0.5,15.5)
ylimits <- c(0,0.3)
samplesize <- 10^5

# simulate a huge sample from Binomial(n, 0.3)
# and plot the probability mass function
Xb <- rbinom(samplesize, 10, 0.3)
hist(Xb, 
     freq = F, 
     col = "DARKRED", 
     xlim=xlimits,
     ylim=ylimits, 
     breaks=thebreaks, 
     xaxt="no",
     xlab="X", 
     main="Probability Mass Function" )
axis(1, at=seq(0,15,1))
axis(1, at=c(11,13,15))

# simulate a huge sample from Poisson(3)
# and plot the probability mass function
Xp <- rpois(samplesize, 3)
par(new=T)
hist(Xp, 
     freq = F, 
     col = "BLUE", 
     xlim=xlimits,
     ylim=ylimits, 
     breaks=thebreaks,
     xaxt="no",
     xlab="X", 
     main="", 
     density=50 )


# tweak the y-axis and replot
ylimits <- c(0,1)
hist(Xb, 
     freq = F, 
     col="DARKRED", 
     xlim=xlimits,
     ylim=ylimits, 
     breaks=thebreaks, 
     xaxt="no",
     xlab="X", 
     main="" )
axis(1, at=seq(0,15,1))
axis(1, at=c(11,13,15))


# now add the cumulative mass function 
# for the Binomial(10, 0.3)
Xbcdf <- ecdf(Xb)
par(new=T)
plot(Xbcdf, 
     xlim=xlimits, 
     ylim=ylimits, 
     main="PDF and CDF (pmf and cmf)", 
     col="DARKRED", 
     xlab="",
     ylab="", 
     cex=1.4)


# and for the Poisson(3)
par(new=T)
hist(Xp, 
     freq = F, 
     col = "BLUE", 
     xlim=xlimits,
     ylim=ylimits, 
     breaks=thebreaks,
     xaxt="no",
     xlab="X", 
     main="", 
     density=50 )
Xpcdf <- ecdf(Xp)
par(new=T)
plot(Xpcdf, 
     xlim=xlimits, 
     ylim=ylimits, 
     main="PDF and CDF (pmf and cmf)", 
     col="BLUE", 
     xlab="",
     ylab="", 
     cex=1.4 )