---
title: "Maximum Likelihood and Maximum A-Posteriori"
output:
  html_document:
    df_print: paged
editor_options:
  chunk_output_type: console
---

```{r include = F, message = F}
library(tidyverse)
```

```{r}
NLL <- function(x, mu, sigma){
  # How consistent are the observations in x
  # with x ~ N(mu, sigma)
  # Here we are computing the Negative Log-Likelihood (NLL)
  # lower NLL => higher likelihood => data is more consistent
  # with model
  -sum(log(dnorm(x, mu, sigma)))
}
```

This function itself returns a function which takes in `params` and returns the NLL of the data `x` (which is fixed once given to `make.my.nll`) under the parameters `params`.

```{r}
make.my.nll <- function(x){
  function(params){
    NLL(x, params[1], params[2])
  }
}

nll <- make.my.nll(x) #deferred evaluation
nll <- function(params){
    NLL(x, params[1], params[2])
  }
```


Let's now generate some data:

```{r}
x <- rnorm(n = 5000, mean = 1, sd = 2)
```

Now, let's find the parameters with which the data is the most consistent:
```{r}
my.nll <- make.my.nll(x)
nlminb(c(0, 1), my.nll)$par
```

This works, but doesn't work great, because `dnorm` gives an imprecise probability density for very small densities. The following uses the `log = TRUE` argument to `dnorm` to get the log-probability.
```{r}
NLL.c <- function(x, mu, sigma){
  # Like NLL, but computing log(dnorm)
  -sum(dnorm(x, mu, sigma, log = TRUE))
}
make.my.nll.c <- function(x){
  my.nll <- function(params){
    NLL.c(x, params[1], params[2])
  }
}

x <- rnorm(n = 1000, mean = 1, sd = 2)
my.nll <- make.my.nll.c(x)
nlminb(c(0, 1), my.nll)$par
```


Now, let's optimize $P(\mu)P(x | \mu, \sigma)$, with the prior $\mu \sim N(0, 1)$

```{r}
make.my.posterior <- function(x, mean.mu, sigma.mu){
  function(params){
    NLL.c(x, params[1], params[2]) - dnorm(params[1], mean = mean.mu, sd = sigma.mu, log = TRUE)
  }
}
```

```{r}
x <- rnorm(n = 1000, mean = 1, sd = 2)
my.post <- make.my.posterior(x, -10, 1)
nlminb(c(0, 1), my.post)$par
```


### Prior-Data conflicts

First, let's make a function that plots the data likelihood (for $\sigma = 1$) as well as the posterior, for different value of $\mu$.

```{r}
plot.lik.posterior <- function(x, mean.mu, sigma.mu){
  possible.mus <- seq(-10, 10, 0.01)
  sigmas <- rep(1, length(possible.mus))
  possible.params <- mapply(c, possible.mus, sigmas, SIMPLIFY = FALSE)
  my.log.post <- make.my.posterior(x, mean.mu, sigma.mu)
  my.nll <- make.my.nll.c(x)
  
  df <- data.frame(mu = possible.mus,
                   log.likelihood = -sapply(X = possible.params, FUN = my.nll),
                   log.posterior = -sapply(X = possible.params, FUN = my.log.post))
  
  
  df.long <- df %>% pivot_longer(cols = c(log.likelihood, log.posterior), names_to = "f", values_to = "log.p")
  
  
  point.ests <- df.long %>% group_by(f) %>% summarize(mu = mu[which.max(log.p)])
  
  ggplot(data = df.long, mapping = aes(color = f)) + 
    geom_line(mapping = aes(x = mu, y = log.p)) + 
    geom_vline(data = point.ests, mapping = aes(xintercept = mu, color = f), linetype = "dashed")
    
    
  
}

```


On one extreme, we have very little data, and a strong prior (a strong prior means that we're pretty sure about the value of e.g. $\mu$).:

```{r}
x <- rnorm(n = 2, mean = 1, sd = 2)
plot.lik.posterior(x = x, mean.mu = 10, sigma.mu = 0.1)
```

On another extreme, we have a lot of data, and a (relatively) weak prior:

```{r}
x <- rnorm(n = 100, mean = 1, sd = 2)
plot.lik.posterior(x = x, mean.mu = 10, sigma.mu = 0.1)
```


An intermediate situation is when we have prior and data that are in conflict that is not resolved. This won't happen with Gaussian priors and likelihoods, but will sometimes happen with other models.

```{r}

lik.plot <- function(probs, lik){
  df <- data.frame(prob = probs,
                   lik = lik)
  ggplot(df, mapping = aes(x = prob, y = lik)) + 
    geom_line() + 
    geom_vline(mapping = aes(xintercept = probs[which.max(lik)]), color = "green")
}

probs <- seq(0.001, 1, 0.002)
prior <- dnorm(x = probs, mean = 0.5, sd = 0.05) + 0.00000001
prior <- prior/sum(prior)

posterior <- function(r, prob, prior){
  lik <- exp(log(probs)*sum(r) + log(1-probs)*(sum(1-r)))
  post <- lik * prior
  post/sum(post)
}

set.seed(100)
r <- rbinom(n = 50, size = 1, prob = 0.95) # mean(r) == 0.98
prior <- dnorm(x = probs, mean = 0.5, sd = 0.05) + 0.0000001
prior <- prior/sum(prior)

lik.plot(probs, posterior(r, probs, prior))

```

