Punto 1

Consideriamo i dati contenuti nel dataframe mcycle all’interno del package MASS e creiamo un primo grafico per visualizzarli. inoltre specifichiamo un seed in modo tale da garantire la ripetibilità dei risultati.

library(MASS)
data(mcycle)
plot(mcycle$times,mcycle$accel)

set.seed(1010)

Punto 2

Suddividiamo in maniera casuale il dataset in due sottoinsiemi, ripsettivamente train e test con una proporzione di 70:30

index<-sample(133,90)
train<- mcycle[index,]
test<- mcycle[-index,]

Piccola modifica al dataset train per rendere gli indici ascendenti (servirà per la costruzione dei grafici)

train[,3]<-index
train<-train[order(train$V3),]
train<-train[,-3]

Visualizziamo in un grafico i 2 sottoinsiemi

plot(train, col="red",ylab= "acceleration")
points(test, col="blue")
legend("topleft",legend = c("Train", "Test"),pch = 21,col = c("blue","red"))

Punto 3

Adattiamo il dataframe train ad un modello lineare, parabolico e cubico

lin<-lm(accel~times,train)
par<-lm(train$accel~train$times+I(train$times^2),train)
cub<-lm(train$accel~train$times+I(train$times^2)++I(train$times^3),train)

plot(train, col="red",ylab= "acceleration")
abline(lin, col="black")
lines(train$times,predict(par,train),col="green",lwd=3)
lines(train$times,predict(cub,train),col="purple",lwd=3)
legend("bottomright",legend = c("K=1","K=2","K=3"),lty =1,lwd=3,col = c("black","green","purple"))

Punto 4

Implementiamo all’interno del comando lm la funzione poly() che ci consente di fare polinomi di grado superiori con facilità

m1<-lm(accel~poly(times,1),train)
m2<-lm(accel~poly(times,2),train)
m3<-lm(accel~poly(times,3),train)
plot(train, col="red",ylab= "acceleration")
lines(train$times,predict(m1,train),col="black",lwd=3)
lines(train$times,predict(m2,train),col="green",lwd=3)
lines(train$times,predict(m3,train),col="purple",lwd=3)
legend("bottomright",legend = c("K=1","K=2","K=3"),lty =1,lwd=3,col = c("black","green","purple"))

Possiamo notare che il grafico che risulta è identico, ora proviamo ad aumentare il numero di parametri

m5<-lm(accel~poly(times,5),train)
m10<-lm(accel~poly(times,10),train)
m15<-lm(accel~poly(times,15),train)
m20<-lm(accel~poly(times,20),train)
plot(train, col="red",ylab= "acceleration")
lines(train$times,predict(m5,train),col="black",lwd=3)
lines(train$times,predict(m10,train),col="green",lwd=3)
lines(train$times,predict(m15,train),col="purple",lwd=3)
lines(train$times,predict(m20,train),col="orange",lwd=3)
legend("bottomright",legend = c("K=5","K=10","K=15","K=20"),lty =1,lwd=3,col = c("black","green","purple","orange"))

Punto 5

Calcoliamo la misura dell’adattamento sul train per i diversi modelli creati grazie alla funzione poly(), per fare ciò useremo il Mean Square Error (MSE)

MSE<-c()
for (k in 1:70){
  
  m<-lm(accel~poly(times,k,raw=T), train)
  MSE[k]<-(sum((train$accel-predict(m,train))^2))/length(train$accel)
}

plot(c(1:70),MSE, type="l",xlab="Flexibility",ylab="MSE")

Come ci aspettavamo, più aumenta la flessibilità (aumentano i parametri) più il modello si adatta ai dati, ottenendo un MSE decrescente

Punto 6

Valutiamo ora l’adattamento dei modelli sul dataset test

plot(test, col="blue",ylab= "acceleration")
lines(test$times,predict(m5,test),col="black",lwd=3)
lines(test$times,predict(m10,test),col="green",lwd=3)
lines(test$times,predict(m15,test),col="purple",lwd=3)
lines(test$times,predict(m20,test),col="orange",lwd=3)
legend("bottomright",legend = c("K=5","K=10","K=15","K=20"),lty =1,lwd=3,col = c("black","green","purple","orange")) 

MSET<-c()
for (k in 1:70){
  
  m<-lm(accel~poly(times,k,raw=T), train)
  MSET[k]<-(sum((test$accel-predict(m,test))^2))/length(test$accel)
}

plot(c(1:70),MSET, type="l",xlab="Flexibility",ylab="MSE")

Punto 7

Confrontiamo i due MSE, quello che ci aspettiamo è di vedere l’MSE del training set che diminuisce con il crescere dei parametri K, mentre l’MSE del test set inizialmente diminuire per poi aumentare con il crescere di K

plot(c(1:70),MSE, type="l",xlab="Flexibility",ylab="MSE",col="blue",lwd=3)
lines(c(1:70),MSET, col="red", lwd=3)
legend("bottomright",legend = c("Train","Test"),lty =1,lwd=3,col = c("blue","red"))

Possiamo constatare che avviene ciò che ci aspettiamo, anche se per veder aumentare l’MSE del test set bisogna arrivare a K elevati (oltre il 50)