Using R: Method of Moments and ML estimators for Beta Binomial Distribution

In this video we derive the method of moments and maximum likelihood estimates from scratch using R.

Help this channel to remain great! Donating to Patreon or Paypal can do this!

[statisticsmatt]

##
## Using R: Method of Moments and ML estimators
## for Beta Binomial Distribution
##
## Follow up video from
## (1) Mean and Variance Beta Binomial
## (2) Method of Moments Beta Binomial
## (3) Maximum Likelihood Beta Binomial
##
## Suggest Reading Article
## "USING THE BETA-BINOMIAL DISTRIBUTION TO ASSESS
## PERFORMANCE OF A BIOMETRIC IDENTIFICATION DEVICE"
##


## Example from Wikipedia
(y <- 0:12)
(freq <- c(3, 24, 104, 286, 670, 1033, 1343, 1112, 829, 478, 181, 45, 7))
x <- rep(y, freq)
(N <- sum(freq))
(n <- max(y))

## Method of Moments
(m1 <- sum(x)/N)
(m2 <- sum(x^2)/N)
(a <-
(b <-

bb <- function(x, n, a, b) choose(n, x) * beta(x + a, n - x + b) / beta(a, b)

tmp<-barplot(table(x), main="Freq overlayed with Beta Binomial Estimate")
lines(tmp, N*bb(0:12, 12, a, b), type='b')



### MLEs ##
## partial derivatives. See video 3 above.
ll.a <- function(a, b) + sum(digamma(a+x))
ll.aa <- function(a, b) + sum(trigamma(a+x))
ll.b <- function(a, b) + sum(digamma(b+n-x))
ll.bb <- function(a, b) + sum(trigamma(b+n-x))
ll.ab <- function(a, b)

## Use MoM as initial estimates
a2<-a;b2<-b
## loop for the iterative process
for(i in 1:50) {
parm<-matrix(c(a2, b2), nrow=2)
f <- matrix(c(ll.a(a2, b2), ll.b(a2, b2)), nrow=2)
J <- solve(matrix(c(ll.aa(a2, b2), ll.ab(a2, b2), ll.ab(a2, b2), ll.bb(a2, b2)), nrow=2, ncol=2))
f2 <- parm - J%*%f
a2<-f2[1, ];b2<-f2[2, ]
}

c("MM"=a, "MLE"=a2, "MM"=b, "MLE"=b2)

## estimates using each type of estimator
round(rbind("observed"=freq, "m.m."=N*bb(0:12, 12, a, b), "mle"=N*bb(0:12, 12, a2, b2), "Bin"=N*dbinom(0:12, 12, a/(a+b))), 1)

tmp<-barplot(table(x), main="Freq overlayed with Binomial and Beta Binomial Estimates", ylim=c(0, 1400))
lines(tmp, N*bb(0:12, 12, a, b), type='b', lwd=2)
lines(tmp, N*dbinom(0:12, 12, a/(a+b)), type='b', col='red', lwd=2)
legend("topright", c('Beta Binomial', 'Binomial'), lty=1, col=c('black', 'red'))


#### Maximized Log Likelihood ####
## Beta Binomial
(LL <- sum(log(bb(x, n, a2, b2))))
## Binomial
(p.hat <- a2/(a2+b2))
(LL2 <- sum(log(dbinom(x, n, p.hat))))

#### AIC (2k - 2 log Likelihood) ####
c("Beta Binomial"=2*2 - 2*LL, "Binomial"=2*1 - 2*LL2)

[kirillivanov]

Great video, could you also show ML parameter estimation for beta distribution using a much simpler "grid search" approach? Thanks

[kirillivanov]

Could you please also do a video of estimating beta distribution parameters via "matching quantiles" or "quantile estimator" method? I'm pretty sure we should have a closed-form solution (or at least be able to solve a non-linear eq.) because we have 2 unknowns and 2 equations. Thanks in advance