Reputation: 47
I need to find the maximum likelihood estimate for a vector of binomial data.
one like this:
binvec <- rbinom(1000, 1, 0.5)
I tried to first create the function and then optimize it with optim()
.
llbinom <- function(theta, x) {return(sum(dbinom(x=x,
size = theta[1],
prob = theta[2],log=TRUE)))}
mybin <- optim(c(0.5,0.5),fn=llbinom,x=binvec)
mybin
I do get some result but also error messages that NaN
s are being produced and the function cannot be evaluated at the initial parameters. I constructed it from an example that works for normally distributed data and believe I made a mistake in converting it.
Here, the original code I got:
ll <- function(theta,x) {return(-sum(dnorm(x=x,
mean=theta[1],sd=theta[2],log=TRUE)))}
mle <- optim(c(5,3),fn=ll,x=binvec)
Upvotes: 1
Views: 1282
Reputation: 226077
Several issues here.
optim()
minimizes by default unless you set the control parameter fnscale=-1
, so you need to define a negative log-likelihood function)size
parameter must be an integersize
parameter from data (this is often done using N-mixture models, if you want to read up on the technique: e.g. see the unmarked
package); usually the number of trials is assumed to be known. So I would tryllbinom <- function(theta, x) {return(
-sum(dbinom(x=x,
size = 1,
prob = theta[1],log=TRUE)))}
mybin <- optim(c(0.5),fn=llbinom,x=binvec)
There are plenty of reasons to do this the Hard Way (numerically); if you really only need to find the MLE of the probability of a single binomial sample x
(independent observations with the same probability of success out of s
trials), the analytical solution is sum(x)/sum(s)
...
Upvotes: 2