```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, fig.align="center")
options(width=110)
```

## Maximum Likelihood Estimation

#### Statistics for Laboratory Scientists ( 140.615 )

### Example from class: Binomial

We have 22 successes out of 100 trials. Some Binomial probabilities.

```{r}
dbinom(22,100,prob=22/100)
dbinom(22,100,prob=0.10)
dbinom(22,100,prob=0.20)
dbinom(22,100,prob=0.25)
p <- seq(0,1,by=0.01)
p
plot(p,dbinom(22,100,prob=p))
```

That is the likelihood function!

```{r}
curve(dbinom(22,100,x),from=0,to=1,n=251,ylab="likelihood")
abline(v=22/100,lwd=2,col="blue")     
curve(dbinom(22,100,x),from=0.1,to=0.4,n=251,ylab="likelihood")
abline(v=22/100,lwd=2,col="blue")     
```

The log likelihood function from class.

```{r}
ll <- function(x,n,p){
  x*log(p)+(n-x)*log(1-p)+log(choose(n,x))
}
```

The log likelihood for some values of p.

```{r}
p <- seq(0.1,0.4,by=0.01)
p
ll(x=22,n=100,p=p)
curve(ll(22,100,x),from=0.1,to=0.4,n=251,ylab="log-likelihood")
abline(v=22/100,lwd=2,col="blue") 
```

You can also use the dbinom() function with argument log=TRUE.

```{r}
dbinom(22,100,p,log=TRUE)
curve(dbinom(22,100,x,log=TRUE),from=0.1,to=0.4,n=251,ylab="log-likelihood")
abline(v=22/100,lwd=2,col="blue")     
```

To find the maximum likelihood estimate you do not need the constant in the log likelihood function, i.e., the last part that does not depend on p.

```{r}
ll.2 <- function(x,n,p){
  x*log(p)+(n-x)*log(1-p)
}
curve(ll.2(22,100,x),from=0.1,to=0.4,n=251,ylab="log-likelihood (version 2)")
abline(v=22/100,lwd=2,col="blue") 
```

### Example from class: Poisson log likelihood

We observe the following outcomes from a Poisson distribution.

```{r}
dat <- c(1,1,0,1,2,5,2,4,0,0,1,0,3,1,0,0,1,1,3,4)
dat
table(dat)
length(dat)
mean(dat)
```

The log likelihood function.

```{r}
ll <- function(dat,lambda){
  -length(dat)*lambda+sum(dat)*log(lambda)-sum(log(gamma(dat+1)))
}
```

The log likelihood for some values of lambda.

```{r}
ll(dat,seq(0.5,4.5,by=0.1))
```

The plot.

```{r}
curve(ll(dat,x),from=0.5,to=4,xlab=expression(lambda),ylab="log likelihood",main="n=20, mean=1.5")
abline(v=mean(dat),lwd=2,col="blue")
```

You can also use the dpois() function with argument log=TRUE.

```{r}
dpois(dat,lambda=1.5,log=TRUE)
sum(dpois(dat,lambda=1.5,log=TRUE))
```