Load packages and set plot theme

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\))

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

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\)

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

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\)

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

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\)

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 0.0078146