```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library("knitr")
library("dplyr")
library("ggplot2")
library("flextable")
library("data.table")
```

### Read the data

```{r poly_data}
poly_data <- fread("../data/poly_data.txt")
```

#### Train/test split + model fit

We define a function that splits the data in training/test sets and then fits the polynomial regression model.
!! Important !!: we specify the **complexity parameter $\delta$ (d)** (for degree of polynomial).

```{r split, echo=FALSE}
## df = data frame with the data
## n = size of the test set
## d = complexity parameter
cross_validate <- function(df, n, d) {
  
  ## split train/test
  vec <- sample(nrow(df), n)
  test <- df[vec, -1]
  train <- df[!vec, -1]
  
  fit <- lm(y ~ poly(x = x, degree = d, raw = TRUE), data = train)
  
  predictions <- predict(fit, test, interval="none", type = "response", na.action=na.pass)
  preds = cbind(test, predictions)
  
  temp <- preds |>
  mutate(squared_error = (y-predictions)^2) |>
  summarise(rmse = sqrt(mean(squared_error)), r = round(cor(y, predictions, method = "pearson"),3), pct_error = 100*rmse/abs(mean(train$y)))
  
  return(temp)
}
```

We try this function once, manually:

```{r}
cross_validate(df = poly_data, n = 4, d = 1) |>
  regulartable() |>
  autofit()
```

### Cross-validation

Now we try different values for the complexity parameter, so to **fine-tune** it.
We define a vector of values for $\delta$; in addition, instead of running the model only once for each value of $\delta$,
we do this multiple times (replicates), each time on a different random split of the data.

```{r}
dval = seq(1,9) ## values of the hyperparameter to be tested
nrep = 50 ## n. of replicates for each value of d
ntest = 5 ## size of the test set

res = data.frame(NULL)

for (d in dval) {
  
  print(paste("trying complexity value = ", d))
  for (i in 1:nrep) {
    
    # print(paste("replicate n.", i))
  
    temp = cross_validate(df = poly_data, n = ntest, d = d)
  
    temp$rep = i
    temp$test_size = ntest
    temp$d = d
  
    res = rbind.data.frame(res, temp)
  }
}

```

```{r}
res$d = factor(res$d)
ggplot(res, aes(x = d, y = log(rmse))) + geom_jitter(aes(color=d)) + geom_boxplot(aes(color=d))
```

```{r}
temp <- res |>
  group_by(d) |>
  summarise(avg = mean(rmse))

temp |>
  regulartable() |>
  autofit()
```

```{r}
temp[which.min(temp$avg),]
```
