Binomial Likelihood

Example 1: Sampling Distribution and Likelihood functions In an experiment, 10 seed are planted. 7 are observed to germinate and produce seedlings.

1. What is an appropriate distribution to use to model this process? Why?

Solution: Binomial distribution as it models numbers of successes in n trials.

2. Plot the Sampling distribution assuming 10 seeds are planted and probability of germination is 0.5: Calculate P(x|p) for x=0,1,2,3,4…10 assuming p->0.5 and plot Probability Mass Function

Approximate solution via simulation then exact solution

x<-rbinom(n=10000,size=10,prob=0.5)
xProb<-hist(x,breaks=c(0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5),plot=T,probability=TRUE)

plot(0:10,dbinom(x=0:10,size=10,prob=0.5),ylab="Prob",xlab="Number of Successes",ylim=c(0,0.4))
mtext(side=3,"Sampling Distribution",line=1)

3. Plot the Likelihood function given 7 of the 10 seeds germinated. Estimate the mle for p?

plot(seq(from=0,to=1,length=100),dbinom(x=7,size=10,prob=seq(from=0,to=1,length=100)),
     typ="p",ylab="Likelihood",xlab="P")
mtext(side=3,"Likelihood")

abline(v=0.70)

3. Plot the Likelihood function given 700 of the 1000 seeds germinated. Estimate the mle for p?

plot(seq(from=0,to=1,length=1000),dbinom(x=700,size=1000,prob=seq(from=0,to=1,length=1000)),
     typ="l",ylab="Likelihood",xlab="P")
mtext(side=3,"Likelihood")

abline(v=0.70)

PROBLEM 1: ************************************************************************** Make a figure with 6 panels that show the likelihoods for the following samples (seedlings,seeds): (3,4),(6,8),(12,16),(24,32),(300,400), and (3000,4000). Hint: Use par(mfrow=c(3,2)) to set up a figure with six panels

Solution:

pVect<-seq(0,1,length=1000)
parOld<-par(mfrow=c(3,2),oma=c(1,1,1,1),mar=c(1,1,1,1)+0.5)
plot(pVect,dbinom(x=3,size=4,prob=pVect),typ="l",ylab="Likelihood",xlab="P")
plot(pVect,dbinom(x=6,size=8,prob=pVect),typ="l",ylab="Likelihood",xlab="P")

plot(pVect,dbinom(x=12,size=16,prob=pVect),typ="l",ylab="Likelihood",xlab="P")
plot(pVect,dbinom(x=24,size=32,prob=pVect),typ="l",ylab="Likelihood",xlab="P")

plot(pVect,dbinom(x=300,size=400,prob=pVect),typ="l",ylab="Likelihood",xlab="P")
plot(pVect,dbinom(x=3000,size=4000,prob=pVect),typ="l",ylab="Likelihood",xlab="P")

par(parOld)

Next moving from likelihood to Bayesian. Example 2: A Beta pdf is often used as a prior for a Binomial. Why? Use a Beta(1,1) and Beta(10,10) as a prior and use a grid approximation to calculate the posterior distribution for problem 1.iii and for problem 2.

First, plotting two separate beta priors: A flat prior and a diffuse prior

a<-1; b<-1
bMean<-a/(a+b)
bVar<-a*b/((a+b)^2 * (a+b+1))
  
pVect<-seq(from=0, to=1, length=1000)
betaDensity<-dbeta(pVect, shape1=1,shape2=1)
plot(pVect,betaDensity,type='l')

pVect<-seq(from=0, to=1, length=1000)
betaDensity<-dbeta(pVect, shape1=10,shape2=10)
plot(pVect,betaDensity,type='l')

Now, plotting the posterior distributions by modifying the likelihood code above to give a grid approximation to posterior probability distributions for Beta(1,1) then Beta(10,10) priors:

a<-b<-1 # Beta parms
parOld<-par(mfrow=c(3,2),oma=c(1,1,1,1),mar=c(1,1,1,1)+0.5)
plot(pVect,
     dbinom(x=3,size=4,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=3,size=4,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")
plot(pVect,
     dbinom(x=6,size=8,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=6,size=8,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

plot(pVect,
     dbinom(x=12,size=16,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=12,size=16,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")
plot(pVect,
     dbinom(x=24,size=32,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=24,size=42,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

plot(pVect,
     dbinom(x=300,size=400,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=300,size=400,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")
plot(pVect,
     dbinom(x=3000,size=4000,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=3000,size=4000,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

par(parOld)

Recalculating with second prior…

a<-b<-10 # Beta parms
parOld<-par(mfrow=c(3,2),oma=c(1,1,1,1),mar=c(1,1,1,1)+0.5)
plot(pVect,
     dbinom(x=3,size=4,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=3,size=4,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")
plot(pVect,
     dbinom(x=6,size=8,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=6,size=8,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

plot(pVect,
     dbinom(x=12,size=16,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=12,size=16,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")
plot(pVect,
     dbinom(x=24,size=32,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=24,size=42,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

plot(pVect,
     dbinom(x=300,size=400,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=300,size=400,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")
plot(pVect,
     dbinom(x=3000,size=4000,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=3000,size=4000,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

par(parOld)

PROBLEM 2: ************************************************************************** Make a figure with 6 panels that show the prior, likelihood, and posterior for the following samples (seedlings,seeds):(1,10),(10,100),(100,1000). Choose both a flat prior (e.g., Beta(1,1) ) and a strong prior using the Beta distribution.

Solution

a<-10; b<-10 # FLAT Beta parms
seeds<-1000; seedlings<-100
pVect<-seq(0,1,length=1000)

parOld<-par(mfrow=c(3,1),oma=c(1,1,1,1),mar=c(1,1,1,1)+0.5)
plot(pVect,dbinom(x=seedlings,size=seeds,prob=pVect),typ="l",ylab="Likelihood",xlab="P")
plot(pVect,dbeta(pVect,shape1=a,shape2=b),typ="l",ylab="Likelihood",xlab="P")

plot(pVect,dbinom(x=seedlings,size=seeds,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)/sum(dbinom(x=seedlings,size=seeds,prob=pVect)*dbeta(pVect,shape1=a,shape2=b)),
     typ="l",ylab="Likelihood",xlab="P")

par(parOld)