Prova 1 - Aprendizagem de Máquina

Author

Joana D’arc Nunes da Silva, Matrícula: 20180078535

Published

14 Aug, 2024

Questão 1:

Disserte sobre o processo de treinamento de modelo de regressão em aprendizagem de máquina. Explique cada um dos passos considerando a imagem que segue:
Não esqueça de explicar:
1. O que é a validação cruzada e qual a diferença entre o leave-one-out e o k-fold cross-validation;
2. O que é o risco preditivo do modelo e qual seu estimador;
3. A importância de se ter uma base de dados de teste para avaliação final do modelo.

Resposta: No processo de treinamento de um modelo de regressão em aprendizagem de máquina, primeimeiramente fazemos o data spliting dos dados, ou seja, a divisão dos dados que usualmente é dividida em treino, validação e teste, onde no conjunto de treino realizamos o ajuste do modelo que melhor se adaptou aos dados, no conjunto de validação fazemos a tunagem de hiperparâmetros do modelo, e no conjunto de teste é onde avaliamos o desempenho do modelo final escolhido por meio do risco preditivo. A validação cruzada consiste na divisão dos dados em treino e validação, em que no leave-one-out cross validation a cada iteração tiramos de fora uma única observação para teste e treinamos o modelo com as observações que ficaram, assim consequentemente, iremos ter \(n\) modelos ajustados, em que \(n\) é o número de observações total na base de dados. Já no k-fold cross validation iremos ter \(k\) modelos finais ajustados, em que \(k\) é o número de lotes, onde a cada iteração um determinado lote será utilizado para teste e o restante das observações serão usadas para treinar o modelo. O risco preditivo avalia o desempenho do modelo ajustado, no qual o Erro Quadrático Médio (EQM) geralmente é utilizado como medida para avaliar este desempenho, onde o EQM é a média do quadrado dos erros entre os valores observados e os valores preditos pelo modelo. É importante termos uma base de dados somente de teste para avaliação do modelo final, pois é nela que verificamos a capacidade de generalização do modelo.

Questão 2:

Considere o modelo de regressão real definido pela equação abaixo:
\[r(x) = 2.76 + 0.5x_1 - 0.75x_2 + 0.5x_3 - 0.75x_4 + x_5 + \sum_{i=6}^{30} 0x_i + \epsilon,\]
em que \(\epsilon \sim N(0,0.5^2)\) e \(x_i \sim N(0,1)\), \(\forall{i} = 1, ..., 30\). Treine um modelo de regressão linear múltipla e estime o risco preditivo do modelo.

Dicas:
1. Considere uma base de dados com \(n=1000\) observações e \(31\) colunas;
2. Note que \(X\) será uma matriz de features (recursos/covariáveis) com \(31\) colunas, sendo a primeira coluna composta por \(1’s\) (por conta do intercepto) e as demais colunas compostas por valores aleatórios de uma distribuição normal padrão, conforme mencionado anteriormente, i.e., \(x_i \sim N(0,1), \forall{i}\);
3. Note que não há hiperparâmetros a serem ajustados, diferentemente da regressão polinomial que vimos em sala de aula;
4. Não havendo hiperparâmetros, você precisa apenas dividir a base de dados entre treino e teste, isto é, realizar o hold-out;
5. Ajuste o modelo na base de dados considerando o conjunto de treino;
6. Avalie o modelo na base de teste.

Resposta:

Code

rm(list=ls())
library(tidymodels)

── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──

✔ broom        1.0.6     ✔ recipes      1.1.0
✔ dials        1.3.0     ✔ rsample      1.2.1
✔ dplyr        1.1.4     ✔ tibble       3.2.1
✔ ggplot2      3.5.1     ✔ tidyr        1.3.1
✔ infer        1.0.7     ✔ tune         1.2.1
✔ modeldata    1.4.0     ✔ workflows    1.1.4
✔ parsnip      1.2.1     ✔ workflowsets 1.1.0
✔ purrr        1.0.2     ✔ yardstick    1.3.1

── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ purrr::discard() masks scales::discard()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ recipes::step()  masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.

Code

library(rsample)
library(patchwork)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ lubridate 1.9.3     ✔ stringr   1.5.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ readr::col_factor() masks scales::col_factor()
✖ purrr::discard()    masks scales::discard()
✖ dplyr::filter()     masks stats::filter()
✖ stringr::fixed()    masks recipes::fixed()
✖ dplyr::lag()        masks stats::lag()
✖ readr::spec()       masks yardstick::spec()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

tidymodels::tidymodels_prefer()

i = 31 
n = 1000
# Fixando a semente
set.seed(0)

# Matriz com as features
X <- cbind(1, matrix(rnorm(n*(i-1), mean = 0, sd = 1), n, i-1))

# Vetor com os coeficientes
coeficientes <- c(2.76, 0.5, -0.75, 0.5, -0.75, 1.0, rep(0,i-6))

# Modelo de regressão real
y <- X %*% coeficientes + rnorm(n, mean = 0, sd = 0.5)

# Juntando as variáveis em um dataframe
dados <- as.data.frame(X)
names(dados) = paste0("x", 0:30)
dados = dados %>% mutate(y = y)

#  Realizando o Hold-out
dados_split <- rsample::initial_split(dados, prop = 0.8, strata = y)
treino <- rsample::training(dados_split)
teste <- rsample::testing(dados_split)

# Definindo o modelo
modelo = lm(y ~ ., data = treino)
modelo


Call:
lm(formula = y ~ ., data = treino)

Coefficients:
(Intercept)           x0           x1           x2           x3           x4  
   2.760087           NA     0.497702    -0.727672     0.504867    -0.763521  
         x5           x6           x7           x8           x9          x10  
   1.022367     0.001728    -0.025276    -0.007400    -0.005678     0.001841  
        x11          x12          x13          x14          x15          x16  
   0.003786    -0.001997    -0.034929    -0.007591     0.011801    -0.005556  
        x17          x18          x19          x20          x21          x22  
   0.001710    -0.008173     0.005426    -0.002834    -0.036208     0.002494  
        x23          x24          x25          x26          x27          x28  
  -0.010351    -0.006369     0.001925     0.010003     0.015120     0.018765  
        x29          x30  
  -0.027867     0.020564

Code

# Fazendo predições com o modeo
y_pred = predict(modelo, newdata = teste)
y_pred

          1           2           3           4           5           6 
 4.43169184  1.55480077  0.58321394  3.28294487  1.10715076  2.24793883 
          7           8           9          10          11          12 
 2.76239035  2.17420385 -0.64577709  2.18540688  3.61188991  1.83812269 
         13          14          15          16          17          18 
-0.19509426 -0.47310159  4.33392194  3.91184945  4.04664312  4.07350094 
         19          20          21          22          23          24 
 2.67550194  2.74354318  2.30787668  2.99719644  4.06180146  4.98428109 
         25          26          27          28          29          30 
 1.97255985  2.40685876  3.84945818  3.26400783  3.93524810  2.60431019 
         31          32          33          34          35          36 
 1.75148043  2.53267140  1.89679593  0.76704382  2.88529071  3.51956809 
         37          38          39          40          41          42 
 4.61011088 -0.30093909  3.16051355 -0.33414965  3.11392948  2.31029979 
         43          44          45          46          47          48 
 3.40209492  3.07513651  0.40894490  3.09182064  5.66155474  3.77924140 
         49          50          51          52          53          54 
 1.88155743  3.53847195  5.02262908  3.90539479  1.78284006  3.90214413 
         55          56          57          58          59          60 
 2.74952173  2.51894698  2.34064720  2.77557124  3.75745748  3.98906013 
         61          62          63          64          65          66 
 1.40253330  1.91924529 -0.79989135  2.70796960  1.87466877  3.27908931 
         67          68          69          70          71          72 
 3.85669475 -0.18680121  2.08917390  1.58887926  1.75962048  3.98413842 
         73          74          75          76          77          78 
 4.94704622  3.61666501  2.32677637  3.44780669  4.65116422  4.39840237 
         79          80          81          82          83          84 
 5.32659388  3.21367040  3.70371305  2.84711168  2.59405022  0.85258378 
         85          86          87          88          89          90 
 2.22576095  0.86030876  5.68685516  3.69887912  0.61845638  3.82102370 
         91          92          93          94          95          96 
 3.10254720  2.37453010  3.30699886  2.47185871  0.49771366  3.10893134 
         97          98          99         100         101         102 
 3.44974730  3.44880695  1.97892613  0.41877799  3.66262666  2.36338256 
        103         104         105         106         107         108 
 2.66396524  4.16364559  2.77282853  4.70673855  5.64936165  2.67999455 
        109         110         111         112         113         114 
 0.30075417  4.40404128  2.19083999  4.99711531  1.81496583  4.00062727 
        115         116         117         118         119         120 
 0.24504624  3.21695614  5.42929550  5.20149864  1.14399461  2.77121367 
        121         122         123         124         125         126 
 5.37439074  4.42327718  0.33076972  1.53292177  2.98323839  3.02089628 
        127         128         129         130         131         132 
 3.66131307  1.25517683  0.53217321  2.86672206 -0.26265836  3.21004464 
        133         134         135         136         137         138 
 3.76927538 -0.55804381 -0.66735396  3.61459400  3.84630572  0.78137432 
        139         140         141         142         143         144 
 4.79431377  4.72058676  1.90344971  3.05031786  2.12123880  0.30174653 
        145         146         147         148         149         150 
 2.96243379  4.97662687  5.68857674 -0.29595233  0.03358453  4.67535351 
        151         152         153         154         155         156 
 3.68771562  2.28397584  1.45727491  2.30751475  4.49802780  2.87812786 
        157         158         159         160         161         162 
 2.34226610  1.57795637  2.69270573 -0.25542712  2.83455261  5.52645521 
        163         164         165         166         167         168 
 2.57021059 -0.91009938 -0.37417016  5.97987430  2.65388582  0.20987721 
        169         170         171         172         173         174 
 4.47013287  2.52142965  6.25943708  2.71299176  5.24500880  1.09279630 
        175         176         177         178         179         180 
 3.18712819  4.76238700  3.75208269  5.36064677  3.11957833  2.99423386 
        181         182         183         184         185         186 
 2.27937563  1.94335225  1.62838134  3.62098415  0.05105675  0.78744093 
        187         188         189         190         191         192 
 5.07171048  4.43894290  8.17489107  5.15438418  0.20867085  3.39240722 
        193         194         195         196         197         198 
 1.56619839  3.98145540  2.45399200  1.19442962  1.40107252  2.36003400 
        199         200 
-0.61627847  3.92044659

Code

# Erro Quadrático Médio
eqm = mean((teste$y - y_pred)^2)
eqm

[1] 0.2642188

O modelo apresentou um erro quadrático médio relativamente baixo de 0.2642.

Questão 3:

Considere o modelo de regressão real dado por: \[r(x) = 1.6 + 5\sin(x) - 8x^2 + \epsilon,\]
em que \(x \sim U(0,20)\) e \(\epsilon \sim N(0,1)\) (normal padrão). Treine modelos de regressão polinomial com o grau do polinômio \(p=1,2,3\) e estime o risco preditivo de cada um dos modelos.

Dicas:
1. Considere \(n = 10000\) observações;
2. Não é necessário fazer fazer cross-validation;
3. Apenas considere treino e teste;
4. Ajuste cada um dos modelos no treino e estime o risco preditivo do modelo de cada um dos modelos no conjunto de teste.

Interprete o resultado obtido.

Resposta:

Code

rm(list=ls())
# Clarregando as bibliotecas
library(tidymodels)
library(ggplot2)
library(tibble)
library(rsample)
library(patchwork)
library(tidyverse)

tidymodels_prefer()

# Modelo real que não conhecemos
random_real <- function(n){
  x <- runif(n = n, min = 0, max = 20)
  # Essa é a regressão real
  y <- 1.6 + 5*sin(x) - 8*(x^2) + rnorm(n, mean = 0, sd = 1)

  tibble(x = x, y = y)
}

# Funções para treinar
treinar <- function(p, dados_treino){
  # Criando uma regressao polinomial
  lm(data = dados_treino, y ~ poly(x, degree = p, raw = TRUE))
}

# Função para testar
teste <- function(modelo, dados_teste){
  # Calculando o erro quadrático médio
  dados_teste %>%
    mutate(
      y_hat = predict(modelo, newdata = dados_teste) # Aqui calculo y chapeu
    ) %>%
    summarise(
      eqm = mean((y - y_hat)^2)
    )
}

# Função para avaliar o modelo
avaliacao_eqm <- function(p, dados_treino, dados_teste){

  # Treinando o modelo
  ajuste <- treinar(p, dados_treino)

  # Calculando o erro quadrático médio
  eqm <- teste(ajuste, dados_teste)

  list(ajuste = ajuste, eqm = as.numeric(eqm[1,1]))
}

Code

# Fixando semente
set.seed(0)

# Gerando a base total de dados com 10000 observações
dados <- random_real(n = 10000)

#  Realizando o Hold-out
dados_divididos <- rsample::initial_split(dados, prop = 0.8, strata = y)
dados_treino <- rsample::training(dados_divididos)
dados_teste <- rsample::testing(dados_divididos)
grau_maximo <- 3

# Avaliando o erro quadrático médio para cada grau
vec_eqm <-
  purrr::map_dbl(
    .x = 1L:grau_maximo,
    .f = \(p) avaliacao_eqm(p = p, dados_treino = dados_treino, dados_teste = dados_teste)$eqm
  )

Code

# Tabela com os EQM's e seus respectivos graus de polinômio
dados_eqm <- tibble(p = 1L:grau_maximo, eqm = vec_eqm)
dados_eqm

# A tibble: 3 × 2
      p     eqm
  <int>   <dbl>
1     1 56140. 
2     2    12.5
3     3    12.3

Code

# Melhor grau do polinômio que tem o menor EQM
melhor_p <- which.min(vec_eqm)
melhor_p

[1] 3

Code

# Ajustando a regressão verdadeira com grau 1 sobre os dados
ajuste1 <- treinar(p = 1, dados_treino = dados_teste)
ajuste1


Call:
lm(formula = y ~ poly(x, degree = p, raw = TRUE), data = dados_treino)

Coefficients:
                    (Intercept)  poly(x, degree = p, raw = TRUE)  
                          538.6                           -160.6

Code

# Ajustando a regressão verdadeira com grau 2 sobre os dados
ajuste2 <- treinar(p = 2, dados_treino = dados_teste)
ajuste2


Call:
lm(formula = y ~ poly(x, degree = p, raw = TRUE), data = dados_treino)

Coefficients:
                     (Intercept)  poly(x, degree = p, raw = TRUE)1  
                          4.1746                           -0.5012  
poly(x, degree = p, raw = TRUE)2  
                         -7.9806

Code

# Ajustando a regressão verdadeira com grau 3 sobre os dados
ajuste3 <- treinar(p = 3, dados_treino = dados_teste)
ajuste3


Call:
lm(formula = y ~ poly(x, degree = p, raw = TRUE), data = dados_treino)

Coefficients:
                     (Intercept)  poly(x, degree = p, raw = TRUE)1  
                        5.265589                         -1.148003  
poly(x, degree = p, raw = TRUE)2  poly(x, degree = p, raw = TRUE)3  
                       -7.900008                         -0.002684

Observa-se que o melhor o modelo de regressão polinomial é o de grau 3, pois apresentou menor erro quadrático médio.

Questão 4:

Considere a base de dados referente à vendas de sorvetes. A base de dados contém as seguintes variáveis:
1. Temperatura: temperatura média do dia;
2. Vendas: quantidade de sorvetes vendidos no dia.

Download: Para baixar os dados, acesse o link e clique em “Download”, no canto superior direito.

Estamos interessados em estimar as vendas de sorvetes dado a temperatura. Dessa forma, considere o número de vendas como sendo o label (variável \(y\)) e as temperaturas como sendo as features (variáveis \(x\)).

Considerando a base de dados fornecida, treine um modelo de regressão polinomial com grau \(p = 1,2,3,4,5\) e estime o risco preditivo do modelo selecionado. Além disso, construa um gráfico do modelo selecionado ajustado aos dados.

Dicas:
1. Considere utilizar um esquema de validação cruzada para selecionar o melhor hiperparâmetro;
2. Com o grau de polinômio escolhido, treine o modelo na base de dados de treino e avalie o risco desse modelo no conjunto de teste;
3. Considere uma divisão inicial (hold-out) na proporção \(80\%\) (treino) e \(20\%\) para teste.

Resposta:

Code

rm(list=ls())
# Carregando os pacotes
library(tidymodels)
library(ggplot2)
library(tibble)
library(rsample)
library(patchwork)
library(tidyverse)
library(corrr)

tidymodels::tidymodels_prefer()

# Carregando a base de dados
dados <- read.csv("C:/Users/gleyc/Downloads/archive (1)/Ice_cream selling data.csv", sep = ",",
                  col.names = c("temperatura","vendas"))

# Renomeando as variáveis
dados <- dados %>% dplyr::rename(x = temperatura, y = vendas)

# Fixando a semente
set.seed(0)

# Realizadon hold-out
dados_split <- rsample::initial_split(dados, prop = 0.8, strata = "y")
treino <- rsample::training(dados_split)
teste <- rsample::testing(dados_split)

# Fazendo validação cruzada
cv <- vfold_cv(treino, v = 5)

# Definindo o modelo
modelo <-
  parsnip::linear_reg() %>%
  parsnip::set_engine("lm") %>%
  parsnip::set_mode("regression")


# Definindo receita
receita <-
  recipe(y ~ ., data = treino) %>%
  recipes::step_poly(
      all_numeric_predictors(),
      degree = tune("p"),
      options = list(raw = TRUE)
    )

# Validacação cruzada
cv <- vfold_cv(treino, v = 5)

# Criando o workflow
wf <-
  workflow() %>%
  add_recipe(receita) %>%
  add_model(modelo)

# Extraindo hiperparametros do modelo
parametros <-
  wf %>%
  extract_parameter_set_dials() %>% # Extraindo os hiperparâmetros do modelo
  extract_parameter_dials("p") %>%
  range_set(range = c(1, 5)) %>%
  parameters()

parametros$id <- "p"

# Tunando o modelo
tunagem <- tune_grid(
  wf,
  resamples = cv,
  grid = tibble(p = 1:5), # Parâmetros que desejamos que ele teste
  metrics = metric_set(rmse)
)

# Coletando as métricas e o melhor grau
collect_metrics(tunagem)

# A tibble: 5 × 7
      p .metric .estimator  mean     n std_err .config             
  <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1     1 rmse    standard   12.7      5   1.41  Preprocessor1_Model1
2     2 rmse    standard    3.66     5   0.407 Preprocessor2_Model1
3     3 rmse    standard    4.51     5   1.35  Preprocessor3_Model1
4     4 rmse    standard    4.20     5   1.53  Preprocessor4_Model1
5     5 rmse    standard    5.02     5   2.35  Preprocessor5_Model1

Code

# Visualizando o melhor grau segundo a métrica RMSE
show_best(tunagem, n = 1, metric = "rmse")

# A tibble: 1 × 7
      p .metric .estimator  mean     n std_err .config             
  <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
1     2 rmse    standard    3.66     5   0.407 Preprocessor2_Model1

Code

# Selecionado o melhor "p"
melhor_p <- select_best(tunagem, metric = "rmse")
melhor_p

# A tibble: 1 × 2
      p .config             
  <int> <chr>               
1     2 Preprocessor2_Model1

O melhor grau de polinômio é 2, segundo a métrica RMSE.

Code

# Finalizando
wf <-
  wf %>%
  finalize_workflow(melhor_p)

# Realizando o ajuste final do modelo
ajuste <- last_fit(wf, dados_split, metrics = metric_set(rmse))

# Olhando o desempenho no teste
collect_metrics(ajuste)

# A tibble: 1 × 4
  .metric .estimator .estimate .config             
  <chr>   <chr>          <dbl> <chr>               
1 rmse    standard        3.08 Preprocessor1_Model1

Code

real_vs_estimado <- collect_predictions(ajuste)

real_vs_estimado %>%
  ggplot(aes(x = .pred, y = y)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "tomato") +
  labs(
    title = "Real vs Estimado",
    subtitle = "Base de teste",
    x = "Previsões",
    y = "Observações"
  )

Code

# Prevendo
dados_qualquer <- sample_n(dados, 10) # Seleciona 10 linhas quaisquer de dados.
wf_modelo <- extract_workflow(ajuste)
valores_previstos <- predict(wf_modelo, new_data = dados_qualquer)

# Visualizando lado a lado, y e valores previstos
dplyr::bind_cols(dados_qualquer, valores_previstos)

             x          y     .pred
1  -2.28826400 18.1239912 14.495141
2   3.70405744 17.8439565 24.299832
3  -2.65228679 20.2796792 18.075166
4   2.07510060  8.1707349  8.807715
5  -0.03389529  0.8976032  2.976671
6   1.24071162  1.2923608  4.581507
7  -1.17312327  6.6891226  6.504845
8   4.13086796 34.5307427 29.942415
9   2.31859124  7.4120940 10.514634
10  3.33593241 26.1047404 19.961191

Code

# Prevendo o número de vendas de sorvetes em uma temperatura de 3 graus celsius
predict(wf_modelo, new_data = tibble(x = 3))

# A tibble: 1 × 1
  .pred
  <dbl>
1  16.4

Estima-se que o número de vendas de sorvetes em uma temperatura de 3 graus celsius é de aproximadamente 16 unidades.

Questão 5:

Considere o modelo de regressão real dado por: \[r(x) = 45 \times tanh\left(\frac{x}{1.9} - 7\right) + 57 + \epsilon,\]
em que \(x\) são observações de uma variável aleatória \(X \sim U(0,18)\) e \(\epsilon \sim N(0,4)\). Considerando um conjunto de dados de 10 mil observações, treine um modelo de regressão polinomial com grau \(p=1, ..., 15\). Estime o risco preditivo do melhor modelo. Construa um gráfico do melhor modelo ajustando aos dados de teste, i.e., \(y\) versus \(\hat{y}\) do conjunto de teste.

Resposta:

Code

rm(list=ls())
# Carregando as bibliotecas
library(tidymodels)
library(ggplot2)
library(tibble)
library(rsample)
library(patchwork)

# Modelo real que não conhecemos
random_real <- function(n){
  x <- runif(n = n, min = 0, max = 18)
  # Regressão real
  y <- 45 * tanh(x/1.9 - 7) + 57  + rnorm(n, 0, 2)

  tibble(x = x, y = y)
}

# Função para treinar o modelo
treinar <- function(p, dados_treino){
  # Criando uma regressao polinomial
  lm(data = dados_treino, y ~ poly(x, degree = p, raw = TRUE))
}

# Função para testar o modelo
teste <- function(modelo, dados_teste){
  # Calculando o erro quadrático médio
  dados_teste %>%
    mutate(
      y_hat = predict(modelo, newdata = dados_teste) # Aqui calculo y chapeu
    ) %>%
    summarise(
      eqm = mean((y - y_hat)^2)
    )
}

# Função para avaliar o erro quadrático médio
avaliacao_eqm <- function(p, dados_treino, dados_teste){

  # Treinando o modelo
  ajuste <- treinar(p, dados_treino)

  # Calculando o erro quadrático médio
  eqm <- teste(ajuste, dados_teste)

  list(ajuste = ajuste, eqm = as.numeric(eqm[1,1]))
}

Code

# Fixando semente
set.seed(0)

# Gerando a base total de dados com 10000 observações
dados <- random_real(n = 10000)

# Realizando o Hold-out
dados_divididos <- rsample::initial_split(dados, prop = 0.8, strata = y)
dados_treino <- rsample::training(dados_divididos)
dados_teste <- rsample::testing(dados_divididos)
grau_maximo <- 15

# Avaliando o erro quadrático médio para cada grau de polinômio
valores <-
  purrr::map_dbl(
    .x = 1L:grau_maximo,
    .f = \(p) avaliacao_eqm(p = p, dados_treino = dados_treino, dados_teste = dados_teste)$eqm
  )

# Melhor grau de polinômio
melhor_p <- which.min(valores)
melhor_p

[1] 14

Code

# Ajustando a regressão verdadeira com melhor_p sobre os dados
ajuste <- treinar(p = melhor_p, dados_treino = dados_teste)

# Visualizando os dados e a regressão real
p <- ggplot(data = dados_teste, aes(x = x, y = y)) +
  geom_point() +
  theme_minimal() +
  ggtitle("Dados reais") +
  stat_function(fun = function(x) 45 * tanh(x/1.9 - 7) + 57, col = "red", linewidth = 1.2)

# Visualizando o ajuste
p_ajuste <-
  p +
  geom_function(fun = function(x) predict(ajuste, newdata = tibble(x = x)), col = "blue", linewidth = 1.2) +
  ggtitle("Ajuste") +
  labs(x = "x", y = "y")

# Visualizando os EQM's
p_eqm <-
  tibble(p = 1L:grau_maximo, eqm = valores) %>%
  ggplot(aes(x = p, y = eqm)) +
  geom_line(linetype = "dotted", linewidth = 1) +
  annotate("point", x = melhor_p, y = valores[melhor_p], col = "red", size = 4, alpha = 0.7) +
  theme_minimal() +
  ggtitle("EQM") +
  labs(x = "p", y = "EQM") +
  scale_y_log10()

p + p_ajuste + p_eqm

Nota-se que o melhor modelo de regressão polinomial é o de grau 14. O gráfico acima mostra a regressão real (em vermelho) e o ajuste do modelo de regressão polinomial de grau 14 (em azul) sobre os dados de teste. Além disso, o gráfico de EQM mostra que o modelo de grau 14 é o que apresenta menor Erro Quadrático Médio.