M o d e l l i n o n l i n e a r i e c r o s s v a l i d a ...taufer/Labs/L3-Temperature_max.pdf ·...

21/10/2017 Modelli non lineari e cross validazione (24)

file:///C:/Users/emanuele.taufer/Google%20Drive/2%20CORSI/3%20SQG/Labs/L3-Temperature_max.html#(24) 1/28

Modelli non lineari e cross validazione

Strumenti quantitativi per la gestioneEmanuele Taufer



Metodi per affrontare problemi non-lineari

Regressione polinomiale

Regressione k-NN

Regressione con funzioni costanti (step functions)

Regressione polinomiale locale (regression splines)

Smoothing splines

Regressione locale

Modelli additivi generalizzati: approccio che permette di usarediversi tipi di tecniche simultaneamente



In questa esercitazione ci concentreremo su:

regressione polinomiale

regressione k-NN

in un problema non-lineare di previsione dei picchi nella domanda dienergia.

I due modelli saranno comparati per mezzo di cross validazione



Regressione polinomialeIl metodo più tradizionale per estendere il modello di regressione acontesti non-lineari è quello di usare una funzione polinomiale:

può essere stimato semplicemente con il metodo dei minimiquadrati: infatti è un modello di regressione standard conpredittori

difficilmente si usa poichè il polinomio tende a diventaretroppo flessibile

la stessa tecnica può essere estesa al caso della regressionelogistica

= + + + ⋯ + +yi β0 β1xi β2x2i βdx

di εi

, , , …xi x2i x3

i xdi

d > 4



Esempio: modellare i picchi didomanda di energia

Per operare efficacemente le aziende elettriche devono essere ingrado di prevedere il picco giornaliero di richiesta di energia

La richiesta di energia è misurata in megawatt per ora

Il picco giornaliero è definito come il massimo di domanda nell’arcodelle 24 ore

Le compagnie elettriche monitorano costantemente il fenomenoattraverso modelli di regressione



Fattori che determinano la domandadi energia

Meteo-sensibili:

Il fattore principale in questo caso è la temperatura (collegata all’usodi impianti di riscaldamento e di condizionamento)

Non meteo-sensibili

Il fattore principale in questo caso è il giorno della settimana: feriale,sabato o domenica (collegato all’uso di frigoriferi, luci, computer,impianti industriali)



I dati Temp-2014-F.datOsservazioni giornaliere di una compagnia elettrica nel periodo 1/1/14- 31/12/14

: Picco di domanda in megawatt (Peak)

: Temperatura in F° al momento del picco (Tmax)

: indicatori se sabato o domenica

Y

X1

,X2 X3



Dataset TempTemp<-read.table("http://www.cs.unitn.it/~taufer/Data/Temp-2014-F.dat",header=T,sep="") head(Temp)

Date Tmax Peak Sat Sun 1 2014-01-01 39.74 3928.294 0 0 2 2014-01-02 34.52 4419.036 0 0 3 2014-01-03 40.10 4114.751 0 0 4 2014-01-04 38.84 4096.177 1 0 5 2014-01-05 42.08 3676.929 0 1 6 2014-01-06 47.66 3410.767 0 0

Per comodità, costruiamo la variabile factor Day

Temp$Day<-Temp$Sat+2*Temp$Sun Temp$Day<-factor(Temp$Day, c(0,1,2),c("Wday","Sat","Sun")) str(Temp$Day)

Factor w/ 3 levels "Wday","Sat","Sun": 1 1 1 2 3 1 1 1 1 1 ...



Plot Picco-Temperatura massimaplot(Temp$Tmax,Temp$Peak,pch=3, col="red",lwd=2.5)



Regressione quadratica - output RPer provare a catturare l’andamento non lineare del fenomeno,utilizziamo un modello di regressione lineare con un terminequadratico. Si noti l’uso del wrapper I() nella formula.

fit2<-lm(Peak~Tmax+I(Tmax^2)+Day,data=Temp) summary(fit2)

Call: lm(formula = Peak ~ Tmax + I(Tmax^2) + Day, data = Temp) Residuals: Min 1Q Median 3Q Max -504.76 -115.79 0.63 113.86 695.00 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.123e+04 1.847e+02 60.801 < 2e-16 *** Tmax -2.505e+02 6.032e+00 -41.533 < 2e-16 *** I(Tmax^2) 1.853e+00 4.755e-02 38.968 < 2e-16 *** DaySat -1.207e+02 2.684e+01 -4.496 9.37e-06 *** DaySun -2.624e+02 2.684e+01 -9.775 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 176.7 on 360 degrees of freedom Multiple R-squared: 0.8645, Adjusted R-squared: 0.863 F-statistic: 574.2 on 4 and 360 DF, p-value: < 2.2e-16



Peak osservato e previstoplot(Temp$Tmax,Temp$Peak,pch=3, col="red",lwd=2.5) points(Temp$Tmax,I(predict(fit2,Temp)),col="blue",pch=8,lwd=2)

I valori previsti dal modello fit2 sono in blu.

Si noti che vi sono tre diverse curve corripondenti ai 3 livelli dellavariabile Day.



Validazione dei modelli di regressione.Per la scelta e la comparazione di modelli di regressione è possibileutilizzare delle statistiche tradizionali che tengono conto del numero dipredittori utilizzati:

-corretto

Se i modelli di regressione da comparare sono innestati fra loro, èpossibile ricorrere all’uso dei test F parziali (vedi lezioni precedenti).

R2

RSE



Cross-validazione (CV)Nel caso in cui si voglia comparare un modello di regressione linearecon un modello di diversa natura, ad esempio un modello diregressione non-parametrica, è opportuno ricorrere alla CV.

La CV può, in generale, essere applicata a tutti i modelli di statisticallearning supervisionato e fornisce un criterio unico per il loroconfronto.

In questa esercitazione useremo

LOOCV

k-fold CV

Nota: poiché in questa esercitazione si userà anche la regressione k-NN, per evitare confusione, parleremo d’ora in poi di m-fold CV



Cross-validazione LOOCV per fit2Il codice seguente calcola una LOOCV utilizzando la funzione glm() ela funzione accessoria cv.glm().

La funzione glm() (generalized lm) può essere utilizzata per stimareun modello di regressione lineare.

library(boot) glm.fit=glm(Peak~Tmax+I(Tmax^2)+Day,data=Temp) cv.err =cv.glm(Temp ,glm.fit)

La stima del test MSE da LOOCV è contenuta nell’oggetto cv.err. Peraccedervi, digitare

cv.err$delta

## [1] 31726.47 31725.16

Dove il primo numero è il valore di test MSE stimato ed il secondouna sua versione corretta per la distorsione (in questo caso quasicoincidono). Vedremo che potranno essere diversi nella m-fold CV.



CV m-fold per fit2Utilizziamo ancora le funzioni glm e cv.glm per ottenere un CV m-foldcon .

cv.err.m =cv.glm(Temp ,glm.fit,K=10)

Da cui

cv.err.m$delta

## [1] 31636.99 31592.38

La stima è molto simile a quella ottenuta con LOOCV

m = 10



Struttura del fenomeno



Regressione k-NNLa regressione quadratica è vincolata alla struttura simmetrica dellaparabola.

Maggior flessibilità potrebbe essere introdotta dall’uso di un terminecubico.

fit3<-lm(Peak ~ Tmax + I(Tmax^2)+I(Tmax^3)+Day,data=Temp)

In quest’esempio utilizziamo una regressione non-parametrica, la k-NN, per cercare di seguire più fedelmente la particolare struttura delfenomeno.



Regressione k-NNPer adattare una regressione k-NN ai dati è necessario costruire unafunzione ad hoc.

La funzione knn.reg() disponibile nel file KNNR.r (scaricabile dallapagina web del corso) produce le previsioni per il vettore di datix.test dato l’input x.train e l’output y.train.

E’ possibile specificare una lista (o anche solo uno) di valori di K daconsiderare

Per indicare ad R dove trovare la funzione knn.reg(), copiare il fileKNNR.r nella directory di lavoro di R e richiamarlo con la funzionesource()

knn.reg <- function(klist,x.train,y.train,x.test) # Regressione k-nearest neighbors # # klist è la lista dei valori K da usare # x.train, y.train: il training set (indipendente-dipendente) # x.test: il test set (o ancora x.train se interessati all'interpolazione) # Output: una matrice di valori previsti per il test set (una colonna per ogni K in klist)

source("KNNR.r")



Input nella funzione knn.reg()In questo caso, la funzione knn.reg(), ci chiede di fornire come inputi dati separati in variabile dipendente -indipendenti, training e, se adisposizione, test.

x.train<-Temp[,c(2,4,5)] y.train<-Temp$Peak

Si noti che per formare x.train non si seleziona la variabile Day, cheha natura factor ed è appropriata per l’utilizzo con lm(), ma siselezionano le due variabili numeriche Sat e Sun.

In questo caso non abbiamo un set di dati test, pertanto utilizziamo lacross-validazione per:

scegliere il valore ottimale di

confrontare la regressione k-NN con quella quadratica

k



CV per la regressione k-NNLa CV può essere fatta utilizzando la funzione knn.reg.cv()disponibile in KNNR.r

knn.reg.cv <- function(klist,x.train,y.train,nfolds)

klist è la lista dei valori K da usare

x.train, y.train: il training set (indipendente-dipendente)

nfolds: il numero di folders per la m-fold CV, indicare n(numerosità dei dati) per ottenere una LOOCV

Output: una matrice di valori previsti per la CV (una colonna per ogniK in klist)



Seleziona il valore k ottimale con 10-fold CV

Il codice seguente, per i valori esegue una 10-foldCV e salva i risultati (ossia le previsioni ottenute dalla CV) nell’oggetto(matrice) p.cv.

I test MSE calcolati successivamente sono nel vettore m.vec.

# Consideriamo valori di k da 1 a 30 klist<-seq(30) ## Ottiene le previsioni da CV (una matrice di dati) p.cv<-knn.reg.cv(klist,x.train,y.train,10) ## Calcola i valori di test MSE stimato per k =1, 2, ... 30 p.cv<-(p.cv-y.train)^2 m.vec<-apply(p.cv, 2, mean)

k = 1, 2, … , 30



Grafico dei test MSE stimati con 10-fold CV

par(mfrow=c(1,1)) plot(klist,m.vec,pch=3, col="red",lwd=2.5)



Determinare il valore minimoUtilizziamo la funzione which.min() per determinare l’elemento(posizione) del vettore m.vec che contiene il valore minimo di MSE

Il modello preferito per la regressione k-NN è quello con k uguale a

which.min(m.vec)

## [1] 9

Il test MSE stimato dalla 10-fold CV è

m.vec[which.min(m.vec)]

## [1] 38040.7



Valori previstiOttenuto il valore ottimale di k, lo utilizziamo per stimare un modelli diregressione k-NN ed ottenere le previsioni.

klist<- which.min(m.vec) fit.knn<- knn.reg(klist,x.train,y.train,x.train)



Valori osservati e valori previsti k-NNpar(mfrow=c(1,1)) plot(Temp$Tmax,Temp$Peak,pch=3, col="red",lwd=2.5) points(Temp$Tmax,fit.knn,col="blue",pch=8,lwd=2)



LOOCV per il modello k-NN sceltoPer confrontare opportunamente fit2 e fit.knn crossvalidiamo laregressione non-parametrica utilizzando LOOCV

klist<- which.min(m.vec) cvr<-knn.reg.cv(klist,x.train,y.train,nrow(Temp)) LOOCV.knn<-mean((cvr-Temp$Peak)^2) LOOCV.knn

## [1] 38552.95



Confronto grafico dei due modellipar(mfrow=c(1,2)) plot(Temp$Tmax,Temp$Peak,pch=3, col="red",lwd=2.5) points(Temp$Tmax,fit.knn,col="blue",pch=8,lwd=2) plot(Temp$Tmax,Temp$Peak,pch=3, col="red",lwd=2.5) points(Temp$Tmax,I(predict(fit2,Temp)),col="blue",pch=8,lwd=2)



Scelta del modelloDalla LOOCV abbiamo una stima del test MSE:

3.1726510^{4} per la regressione quadratica

3.855310^{4} per la regressione k-NN

La regressione quadratica pur non catturando adeguatamente lastruttura del fenomeno sembra sfruttare appieno l’informazione datadal giorno della settimana.

Se si considera anche la facilità di interpretazione, fit2 è decisamentepreferito a fit.kNN.

Un modello più adatto potrebbe essere ottenuto utilizzando le spline(non trattate nel corso)

M o d e l l i n o n l i n e a r i e c r o s s v a l i d a ...taufer/Labs/L3-Temperature_max.pdf ·...

Documents

Transcript of M o d e l l i n o n l i n e a r i e c r o s s v a l i d a ...taufer/Labs/L3-Temperature_max.pdf ·...