# Builds _data/diagnose.RData for scripts/06-diagnose-and-fix.R.
# INSTRUCTOR KEY -- this file is not part of the published site.
#
#  data   n    what is wrong                          fix
#  d1     60   curvature + fan (exponential growth)   log(y)
#  d2     25   nothing (just noisy)                   none
#  d3     70   hump, peak near mean(x); slope ~ 0     y ~ I((x - mean(x))^2)
#  d4     80   strong curvature, x spans 1-1000       log(x)
#  d5     35   no relationship                        y ~ 1
#  d6     50   curvature + mild fan                   y^(1/2)
#  d7     45   annual series: correlated errors;      centre year for the intercept;
#              intercept is the value in year 0       correlation has no fix yet
#  d8     80   nothing -- same model as d4, but       none
#              x only spans 400-900
#  d9     100  curvature + fan, both variables skewed log(y) ~ log(x)
#  d10    70   U shape, low point off-centre (3.5)    y ~ I((x - c)^2), search c
#  d11    40   convex decay + skew                    y^(-1)
#  d12    45   convex decay + fan (exponential decay) log(y)
#
# Pairs worth asking about afterwards:
#   d3 vs d5  -- near-identical summary() output, only the residuals differ
#   d4 vs d8  -- the same model; only the range of x differs

one <- function(seed, n, xgen, ygen) {
  set.seed(seed)
  x <- xgen(n)
  data.frame(x = x, y = ygen(x, n))
}

d1 <- one(102, 60,  function(n) runif(n, 0, 10),
          function(x, n) exp(1 + 0.3*x + rnorm(n, 0, 0.35)))
d2 <- one(101, 25,  function(n) runif(n, 0, 20),
          function(x, n) 4 + 1.5*x + rnorm(n, 0, 8))
d3 <- one(29,  70,  function(n) runif(n, 0, 12),
          function(x, n) 20 - 0.6*(x - 6)^2 + rnorm(n, 0, 2))
d5 <- one(3,   35,  function(n) runif(n, 0, 50),
          function(x, n) 12 + rnorm(n, 0, 3))
d6 <- one(2,   50,  function(n) runif(n, 0, 10),
          function(x, n) (2 + 1.2*x + rnorm(n, 0, 0.8))^2)
d9 <- one(105, 100, function(n) 10^runif(n, 0, 4),
          function(x, n) 2 * x^0.6 * exp(rnorm(n, 0, 0.4)))
d10 <- one(109, 70, function(n) runif(n, 0, 12),
           function(x, n) 5 + 0.4*(x - 3.5)^2 + rnorm(n, 0, 2.5))
d11 <- one(1,  40,  function(n) runif(n, 1, 10),
           function(x, n) 1 / (0.05 + 0.04*x + rnorm(n, 0, 0.015)))
d12 <- one(1,  45,  function(n) runif(n, 0, 10),
           function(x, n) exp(4 - 0.35*x + rnorm(n, 0, 0.3)))

# d4 and d8 share one model and one random stream
set.seed(12)
n  <- 80
xw <- 10^runif(n, 0, 3)
xn <- runif(n, 400, 900)
d4 <- data.frame(x = xw, y = 5 + 30*log(xw) + rnorm(n, 0, 8))
d8 <- data.frame(x = xn, y = 5 + 30*log(xn) + rnorm(n, 0, 8))

# d7: one value per year, in time order, AR(1) errors
set.seed(30)
year <- 1976:2020
e <- as.numeric(arima.sim(list(ar = 0.7), length(year), sd = 3*sqrt(1 - 0.7^2)))
d7 <- data.frame(year = year, y = 30 + 0.4*(year - 1950) + e)

# round for display; keeps every diagnostic the same to two decimals
for (nm in paste0('d', 1:12)) {
  d <- get(nm)
  d$y <- round(d$y, 3)
  if ('x' %in% names(d)) d$x <- round(d$x, 3)
  assign(nm, d)
}

save(list = paste0('d', 1:12), file = here::here('_data', 'diagnose.RData'))