R Estimating parameters of binomial distribution

Question

I'm trying estimate parameters n and p from Binomial Distribution by Maximum Likelihood in R.

I'm using the function optim from stats package, but there is an error.

That is my code:

xi = rbinom(100, 20, 0.5) # Sample
n = length(xi) # Sample size

# Log-Likelihood
lnlike <- function(theta){
log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
(n*theta[1] - sum(xi))*log(1-theta[2])
}

# Optimizing 
optim(theta <- c(10,.3), lnlike, hessian=TRUE)

Error in optim(theta <- c(10, 0.3), lnlike, hessian = TRUE) : function cannot be evaluated at initial parameters

Anyone done this? Which function used?

I don't know how to get to the answer, but I also don't see an error... Could you post it?? Your question should be more focused around how to fix the error over how to get to the solution. (Getting to the solution might be just fixing the error though..) — Cayce K

Ben Bolker Ben Bolker · Accepted Answer · 2016-05-31T18:48:04

tl;dr you're going to get a likelihood of zero (and thus a negative-infinite log-likelihood) if the response variable is greater than the binomial N (which is the theoretical maximum value of the response). In most practical problems, N is taken as known and just the probability is estimated. If you do want to estimate N, you need to (1) constrain it to be >= the largest value in the sample; (2) do something special to optimize over a parameter that must be discrete (this is an advanced/tricky problem).

First part of this answer shows debugging strategies for identifying the problem, second illustrates a strategy for optimizing N and p simultaneously (by brute force over a reasonable range of N).

Setup:

set.seed(101)
n <- 100
xi <- rbinom(n, size=20, prob=0.5) # Sample

Log-likelihood function:

lnlike <- function(theta){
    log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
       (n*theta[1] - sum(xi))*log(1-theta[2])
}

Let's break this down.

theta <- c(10,0.3)  ## starting values
lnlike(c(10,0.3))  ## -Inf

OK, the log-likelihood is -Inf at the starting value. Not surprising that optim() can't work with that.

Let's work through the terms.

log(prod(choose(theta[1],xi))) ## -Inf

OK, we're already in trouble on the first term.

prod(choose(theta[1],xi)) ## 0

The product is zero ... why?

choose(theta[1],xi)
##  [1] 120 210  10   0   0  10 120 210   0   0  45 210   1   0

Lots of zeros. Why? What are the values of xi that are problematic?

## [1]  7  6  9 12 11  9  7  6

Aha! We're OK for 7, 6, 9 ... but in trouble with 12.

badvals <- (choose(theta[1],xi)==0)
all(badvals==(xi>10))  ## TRUE

If you really want to do this, you can do it by brute-force enumeration over reasonable values of n ...

## likelihood function
llik2 <- function(p,n) {
    -sum(dbinom(xi,prob=p,size=n,log=TRUE))
}
## possible N values (15 to 50)
nvec <- max(xi):50
Lvec <- numeric(length(nvec))
for (i in 1:length(nvec)) {
    ## optim() wants method="Brent"/lower/upper for 1-D optimization
    Lvec[i] <- optim(par=0.5,fn=llik2,n=nvec[i],method="Brent",
                     lower=0.001,upper=0.999)$val
}
nvec[which.min(Lvec)]  ## 20
par(las=1,bty="l")
plot(nvec,Lvec,type="b")

R Estimating parameters of binomial distribution

2 Answers