---
title: "When a good fit can be bad"
output:
  html_document:
    highlight: kate
    theme: cosmo
editor_options: 
  chunk_output_type: console
---
Load packages and set plot theme
```{r load-libraries, message=FALSE}
library(tidyverse)
library(ggthemes)

theme_set(theme_few(base_size = 16))

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, cache = TRUE, 
                  message = FALSE, warn = FALSE)
```

Set up the two models ($M_a$ and $M_b$)
```{r setup-models}
m_a <- function(t, a = .5, epsilon = 0) {
  tibble(t = t, y = (1+t)^-a + rnorm(length(t), 0, epsilon))
} 

m_b <- function(t, a = .5, b = .3, c = .7, epsilon = 0) {
  tibble(t = t, y = (b+c*t)^-a + rnorm(length(t), 0, epsilon))
} 
```

Set up the loss functions for each o the two models using RMSE
```{r define-losses}
loss_a <- function(a, observed_data) {
  new_a <- m_a(observed_data$t, a) 
  
  sqrt(sum((new_a$y - observed_data$y)^2) / nrow(observed_data))
}

loss_b <- function(par_list, observed_data) {
  a <- par_list[1]
  b <- par_list[2]
  c <- par_list[3]
  
  new_b <- m_b(observed_data$t, a, b, c) 
  
  sqrt(sum((new_b$y - a_output$y)^2) / nrow(observed_data))
}

```

Simulate samples from Model $M_a$
```{r simulate-a}
a_output <- m_a(1:100, a = .4, epsilon = .1)

ggplot(a_output,aes(x = t, y= y)) + 
  geom_point()
```

Fit model $M_a$ to the simulated data
```{r fit-a}
a_opt <-optim(par = .3, fn = function(x){loss_a(x, a_output)}, 
              method = "Brent", lower = 0, upper = 10)
a_infer <- m_a(t = a_output$t, a = a_opt$par)
```

Plot the estimated model $M_a$
```{r plot-a}
ggplot(a_output,aes(x = t, y= y)) + 
  geom_point() + 
  geom_line(data = a_infer, color = "darkred")
```


Fit model $M_b$ to the simulated data
```{r fit-b}
b_opt <- optim(par = c(.3, .2, .1), fn = function(x){loss_b(x, a_output)})
b_infer <- m_b(t = a_output$t, a = b_opt$par[1], b = b_opt$par[2], c = b_opt$par[3])
```

Plot the estimated model $M_b$
```{r plot-b}
ggplot(a_output,aes(x = t, y= y)) + 
  geom_point() + 
  geom_line(data = a_infer, color = "darkred") + 
  geom_line(data = b_infer, color = "darkblue")
```

The difference in fit between $M_a$ and $M_b$ was `r a_opt$value - b_opt$value`

