# Maximum likelihood for data from a Binomial distribution

I'm trying to follow how maximum likelihood works by using R. I'm following the example here but with some other data.

But I'm confused by the output. How is the MLE so different from what I'd expect which is 0.7 i.e. number of successes / total number of trials?

# MLE for Binomial Distribution
y<-c(0,0,0,1,1,1,1,1,1,1)
n<-length(y)
# formulation for the log likelihood for the binomial
logL <- function(p) sum(log(dbinom(y, n, p)))
# again we can test the function for one value of p
logL(0.8)
#plot logL
p.seq <- seq(0, 0.99, 0.01)
plot(p.seq, sapply(p.seq, logL), type="l")
#optimum:
optimize(logL, lower=0, upper=1, maximum=TRUE)

• I suspect you've gotten "n" wrong; in this case, it looks like you intend your data to be distributed Bernoulli(0.7) (i.e., Binomial with n=1 and p=0.7) and you have 10 observations. If your data really was Binomial with n=10, your 10 observations would be from a Binomial(n=10,p), and since your counts are all 0 or 1, p must be quite low! Try fixing n=1 inside logL and see what you get. Dec 6, 2017 at 17:52
• Thanks! That's what I was after. To clarify then, in this case are my data drawn from multiple Bernoulli distributions? Dec 6, 2017 at 18:08
• Yes. You could look at it as 10 draws from a Bernoulli distribution or 1 draw from a Binomial distribution with n=10; in the latter case, you'd have to sum up the values to get "7", which would be the value of your one observation. If you write your likelihood both ways, you should get identical results from your plot and maximization. Dec 6, 2017 at 18:17
• Could you now answer your own question as you seem to have solved it? so it's not seen as unanswered? Oct 19, 2018 at 22:21

As jbowman noted in the comments, all I need to do is change the value of n from the length of the vector to 1. Then I get the MLE as expected ~ 0.7.

# MLE for Binomial Distribution
y<-c(0,0,0,1,1,1,1,1,1,1)
n<-1 # length(y) here's the change!
# formulation for the log likelihood for the binomial
logL <- function(p) sum(log(dbinom(y, n, p)))
# again we can test the function for one value of p
logL(0.7)
#plot logL
p.seq <- seq(0, 0.99, 0.01)
plot(p.seq, sapply(p.seq, logL), type="l")
#optimum:
optimize(logL, lower=0, upper=1, maximum=TRUE)