Maximum likelihood

class: title-slide, middle

# Maximum likelihood

## Solutions to exercises

.instructors[
  **ECL707/807** - Dominique Gravel & Andrew MacDonald
]

---

# Exercise 2.2. 
## Time to practice : maple distribution at Sutton

![](09_likelihood_solution_files/figure-html/unnamed-chunk-1-1.png)
---

# Exercise 2.2.
## A species distribution problem

**The model**

`$$y \thicksim f(E, \theta)$$`

- f is a probability density function
- `$y$` is the abundance (in stems per quadrat) of sugar maple
- `$E$` is elevation
- `$\theta$` is a set of parameters

**Things to think about :** 
- What are the characteristics of the data ? 
- What is the form of the function ? 
- What is the probability density function ?

---

# Solution 2.2.
## Linear regression

```r
### Likelihood function
ll_fn <- function(a,b,sig,E,obs) {

# Function for the mean
    mu <- a + b*E

# PDF 
    lik <- dnorm(x=obs, mean = mu, sd = sig)
    loglik <- log(lik)

# Return loglikelihood
    sum(loglik)
}
```

---

# Solution 2.2.
## Linear regression

```r
### Try different values 
ll_fn(a=25,b=-0.075,sig=10,E=sutton$y,obs=sutton$acsa)
```

```
## [1] -3275.689
```

```r
ll_fn(a=25,b=-0.05,sig=10,E=sutton$y,obs=sutton$acsa)
```

```
## [1] -2046.312
```

```r
ll_fn(a=30,b=-0.075,sig=10,E=sutton$y,obs=sutton$acsa)
```

```
## [1] -2878.139
```

---

# Solution 2.2.
## Linear regression

```r
### Run a grid search 
a_vec <- seq(15,45,length.out = 100)
b_vec <- seq(-0.1,-0.01,length.out = 100)
res <- matrix(nr=100,nc=100)
for(i in 1:100) {
    for(j in 1:100) {
        res[i,j] <- ll_fn(a=a_vec[i],b=b_vec[j],sig=10,E=sutton$y,obs=sutton$acsa)
    }
}
```
---

# Solution 2.2.
## Linear regression

![](09_likelihood_solution_files/figure-html/unnamed-chunk-5-1.png)

---

# Solution 2.2.
## Logistic regression

```r
### Likelihood function
ll_fn <- function(a,b,sig,E,obs) {

# Function for the mean
    mu <- a + b*E

# logit transformation 
    p = exp(mu)/(1+exp(mu))

# no PDF for logistic regression
    # use the output of the model directly 
    lik <- numeric(length(obs))
    lik[obs==1] = log(p[obs==1])
    lik[obs==0] = log(1-p[obs==0])

# Return loglikelihood
    sum(log(lik))
}
```

---

# Solution 2.2.
## Poisson regression

```r
### Likelihood function
ll_fn <- function(a,b,sig,E,obs) {

# Function for the mean
    mu <- a + b*E

# PDF 
    lik <- dpois(x=obs, lambda = mu)
    loglik <- log(lik)

# Return loglikelihood
    sum(loglik)
}
```