Statistical learning

Statistical learning

Link utili

• mboost

• boosting generale

• plotly

• Grafici

Informazioni utili

• guardo per velocizzare i modelli complessi.

H2O

• plotly (visualizzazione grafica).

Indice:

Lezione 1

Lezione 2

Lezione 3

Lezione 4

Lezione 5

Lezione 6

Lezione 7

Lezione 8

Lezione 9

Lezione 10

Lezione 11

Lezione 12

Lezione 1

Introduzione (con descrizione e codici)

Prima di tutto bisogna inizializzare l’algoritmo (solitamente con la media). Come weak learner oggi prendiamo

−

d d

un albero decisionale e il parametro di complessità è il numero di split e le interazioni sono 1 (vedremo

poi nella seconda lezione il significato di tutto questo). Come weak learner prenderemo alberi con un solo

split. Prendo i residui dei dati con la funzione che ho inizializzato (la media); Scelgo lo split che massimizza

l’eterogeneità tra due gruppi. Ora aggiorno la stima inziale con quella appena trovata (la devo moderare con

ν), analizzo nuovamente i residui e trovo la nuova stima (sempre un albero) ed eseguo l’aggiornamento al

passo successivo. Mi fermo quando non c’è più un cut off che mi fa decrescere l’MSE. Parto quindi da una

stima non informativa e la aggiorno passo a passo.

Con gli alberi decisionali uso ancora il boosting per cercare di migliorare gli errori (attraverso i residui) che

faccio con la stima; Aggiorno ricorsivamente per migliorare gli errori che faccio.

Proviamo a riprodurre l’algoritmo visto con gli alberi.

Gli alberi in una foresta casuale sono identicamente distribuiti mentre, per il boosting, ogni albero prova a

ridurre l’errore facendo l’ensemble degli alberi precedentemente fatti crescere.

Il boosting rischia di andare in overfitting con tanti alberi mentre le random forest no, quindi se trovo il

numero di alberi corretti il boosting è più potente. 1

Il boosting va in overfitting dato che l’MSE di test tende poi a risalire dopo un certo punto. Random Forest

ha una capacità di apprendimento più veloce.

Se non individuiamo il numero di alberi corretti, abbiamo che il boosting va meglio della foresta. Il boosting

ci dice quali e quanti predittori sono stati utilizzati nell’algoritmo.

Problemi di regressione con la perdita di errore quadratico

0

∈

Y X , ..., X

La variabie risposta è e i predittori sono = (X ) ,

R p

1

Il problema di minimizzazione è: {loss[Y,

arg min E f (X)]}

f Y,X

Con la perdita di errore quadratico: 2

−

loss[Y, f f

(X)] = [Y (X)]

|X

f E(Y x).

che ha come soluzione (x) = = ˆ

f f

L’obiettivo è quello di utilizzare il dataset di training per ottenere una stima (x) di (x).

Boosting con la perdita di errore quadratico

ˆ

[0]

f y.

1) Inizializzare la stima di = Lo impostiamo come la media delle y, quando non abbiamo altre

b B >

informazioni. Impostiamo = 0 e fissiamo un numero di interazioni grande con 0.

b

2) Facciamo crescere di 1 e ricaviamo i residui:

ˆ

[b−1]

−

r y f i ..., n

= (x ) per = 1,

i i i

g(x)

3) Fitto il modello ai dati: [b]

→

, x ..., , x ĝ

(r ), (r ) (x)

n n

1 1

4) Aggiorno ˆ ˆ

[b] [b−1] [b]

f f νĝ

(x) = (x) + (x)

< ν <

Dove 0 1 è la step-length factor.

b B B).

5) iteriamo i passi da 2 a 4 finchè non sarà uguale a (b =

g(x).

Torniamo un attimo al passo 3 e osserviamo il modello

g(x) è un modello, un o e deve avere un apprendimento debole, non molto

weak learner base learner

g(x) d

complesso. Possiamo considerare ad esempio come un albero di regressione con splits (profondità

dell’albero). −

d

La complessità del parametro = (d 1) ordine delle interazioni del modello.

→ → d g(x)

albero di regressione con un singolo split = 1 (cosideriamo nelle seguenti analisi = stump.

Stump

Osserviamo un esempio pratico con tanto di codici di ciò che abbiamo appena definito.

Y

= sin(X) +

∼

X U (0, 2π)

∼

N (0, 0.25)

set.seed(123)

n <- 100 #fisso una numerosità

x <- #vettore delle x determinato da un'estrazione casuale da un'uniforme che moltiplica 2

runif(n)*2*pi

y <- + sd=0.5) #definita sopra, da seno più normale di lunghezza pari ad n e con varianz

sin(x) rnorm(n,

#rappresentazione grafica dello scatterplot che definisce le y e le x

plot(y~x) col=2, add=T) #linea che approssima i punti sopra definiti

curve(sin(x), min(x), max(x), 2

2

1

0

y −1

−2 0 1 2 3 4 5 6

x

library(rpart)

#libreria che permette di generare gli alberi di decisione

# step 0

B = 50; nu = 0.1; d = 1 #inizializzo con ciò che mi serve:

#B: numero delle iterazioni

#v: step-length factor ed è un valore compreso tra 0 e 1

#d: parametro di complessità

# step 1

hatf0 <- #calcolo la media delle y

mean(y)

# step 2

r = y - hatf0 #calcolo i residui dati dalla differenza tra veri valori di y e la media di y (ovviamente

# step 3

hatg1 = control=rpart.control(maxdepth = d))

rpart(r~x,

#simple model, spesso definito weak learner o base learner.

#abiamo considerato un albero di regressione con d=1 splits (profondità dell'albero).

#lo definiamo STUMP, cioè un albero di regressione con un solo split d=1.

# step 4

hatfb = nrow=n, ncol=B) #generiamo una matrice composta da soli valori NA con numero di righe

matrix(NA,

hatfb[,1] = hatf0 + nu*predict(hatg1) #consideriamo la prima colonna con valori dati dalla media delle y

#ripeto i procedimenti per b piccolo che va da 2 a B grande, iterando gli step dal passo 2 al passo 4.

3

# step 5

for (b in 2:B){

# step 2

r = y - hatfb[,b-1]

# step 3

hatgb = control=rpart.control(maxdepth = d))

rpart(r~x,

# step 4

hatfb[,b] = hatfb[,b-1] + nu*predict(hatgb)

}

# in questo modo tutte le colonne sono state coperte dai valori grazie al ciclo che ha coperto tutti i v

# estimated f at step bstop

bstop = 40 #stimiamo quindi la funzione f definita sopra considerando un punto b di stop pari a 40. più

plot(y~x) col="red", ylab="y",add=T)

curve(sin(x), min(x),max(x),

ix = #ordinamento in senso crescente la il vettore di valori x per poi rapp

sort(x,index.return=TRUE)$ix

type="s", col="blue")

lines(x[ix],hatfb[ix,bstop],

2

1

0

y −1

−2 0 1 2 3 4 5 6

x

4

Lezione 2

ALS dataset

From Efron and Hastie (2016) Computer-Age Statistical Inference: Algorithms, Evidence and Data Science,

Cambridge University Press, Chapter 17: Random Forests and Boosting.

ALS data represent measurements on patients with amyotrophic lateral sclerosis (Lou Gehrig’s disease). The

goal is to predict the rate of progression of an ALS functional rating score (FRS). There are 1197 training

measurements on 369 predictors and the response, with a corresponding test set of size 625 observations.

Import the data

# import data

als <- Learning/Dataset/ALS.txt", sep="")

read.csv("F:/Statistical

## Warning in read.table(file = file, header = header, sep = sep, quote =

## quote, : line 3 appears to contain embedded nulls

## Warning in read.table(file = file, header = header, sep = sep, quote =

## quote, : line 4 appears to contain embedded nulls

## Warning in scan(file, what = "", sep = sep, quote = quote, nlines = 1,

## quiet = TRUE, : embedded nul(s) found in input

## Warning in scan(file = file, what = what, sep = sep, quote = quote, dec =

## dec, : embedded nul(s) found in input

#1822 osservazioni con 371 predittori

dim(als)

## [1] 630 371

#Split into training and test sets:

# training and test sapendo che la prima colonna ci permette di capire qual è il training e qual è il te

train <- als[als$testset==F,-1]

( n <- )

nrow(train)

## [1] 294

test <- als[als$testset==T,-1]

( m <- )

nrow(test)

## [1] 333

#stimiamo ora il nostro modello:

#Boosting with Squared-Error Loss by using the R package gbm:

# library gbm

library(gbm)

## Loading required package: survival

## Loading required package: lattice

## Loading required package: splines

## Loading required package: parallel

## Loaded gbm 2.1.1

set.seed(123)

Parametri di tuning:

→

d

- = 4 profondità dell’albero. 5

→

ν

- = 0.02 step length factor.

→

B

- = 500 numero di iterazioni.

→

K

- = 10 cross-validazione.

# tuning parameters

d <- 4 #profondità albero

nu <- 0.02 #step length factor

B <- 500 #numero di iterazioni

K <- 10 #cross-validazione (parametro di tuning)

# boosted regression trees fit

fit <- ~ .,

gbm(dFRS

distribution = "gaussian", #di default

data = train,

n.trees = B,

interaction.depth = d,

shrinkage = nu,

bag.fraction = 1, #sto considerando tutte le osservazioni, non faccio un sottocampionamento

cv.folds= K) #inseriamo infine un termine di cross validazione

## Warning in gbm.fit(x, y, offset = offset, distribution = distribution, w =

## w, : variable 33: Uncle..Paternal. has no variation.

print(fit)

## gbm(formula = dFRS ~ ., distribution = "gaussian", data = train,

## n.trees = B, interaction.depth = d, shrinkage = nu, bag.fraction = 1,

## cv.folds = K)

## A gradient boosted model with gaussian loss function.

## 500 iterations were performed.

## The best cross-validation iteration was 1.

## There were 369 predictors of which 3 had non-zero influence.

L’output ci dice che abbiamo 369 predittori dei quali 136 li stiamo utilizzando, dato che hanno influenza non

nulla. Tutto il resto ha influenza nulla.

Test MSE:

# predictions at each step b

#yhats <- sapply(1:B, function(b) predict(fit, newdata = test, n.trees = b))

# genera un vettore di lunghezza che va da 1 a B ed è composto dalla previsione del nostro modello boost

# test MSE at each step b

#MSEs <- apply(yhats, 2, function(yhat) mean( (test$dFRS - yhat)^2))

#plot(1:B, MSEs, xlab="Number of trees", ylab="test MSE", type="l") #graficmente osservo che il boosting

#Quindi vado a vedere il minimo in cosa corrisponde.

#which.min(MSEs)

#il minimo corrisponde a 188 numero di alberi

Early stopping g(X) g(X , X , ..., X X X

Utilizzo uno stump per = ) = , dove è il j-esimo predittore selezionato per lo

∗ ∗

p j j

1 2

split, con j che va da 1 a p. p p

B

X X X X

ˆ ˆ ˆ ˆ ˆ ˆ

[B] [0] [b] [0] [b] [0]

f , ..., X f f , ..., X f νĝ f f

(X ) = + (X ) = + (X ) = + (X )

p p j j j

1 1 j=1 j=1

b=1 b∈B j

6

# bstop

bstop <- method="cv")

gbm.perf(fit,

15000

loss 10000

error

Squared 5000

0 0 100 200 300 400 500

Iteration

#Variable importance (display first 10)

# variable importante

n.trees=bstop)[1:10,]

summary(fit, 7

Family

sd.weight

slope.bp.systolic.slope 0 10 20 30 40 50 60 70

Relative influence

## var rel.inf

## meansquares.slope.fvc.liters meansquares.slope.fvc.liters 72.031225

## Onset.Delta Onset.Delta 24.412255

## last.svc.liters last.svc.liters 3.556519

## Symptom.Speech Symptom.Speech 0.000000

## Symptom.WEAKNESS Symptom.WEAKNESS 0.000000

## Symptom.OTHER Symptom.OTHER 0.000000

## Symptom.Swallowing Symptom.Swallowing 0.000000

## Symptom.GAIT_CHANGES Symptom.GAIT_CHANGES 0.000000

## Symptom.Atrophy Symptom.Atrophy 0.000000

## Symptom.Cramps Symptom.Cramps 0.000000

ˆ

f

Osservo che le funzioni stimate graficamente corrispondono a (X ) sotto rappresentate

j j

# fitted functions for the best 3 predictors

i.var="Onset.Delta", n.trees=bstop)

plot(fit,

## Warning in axis(side = 1, at = 1:499, labels = c("-1000", "-1005",

## "-1007", : larghezza font sconosciuta per il carattere 0x19

## Warning in axis(side = 1, at = 1:499, labels = c("-1000", "-1005",

## "-1007", : larghezza font sconosciuta per il carattere 0x17

## Warning in axis(side = 1, at = 1:499, labels = c("-1000", "-1005",

## "-1007", : larghezza font sconosciuta per il carattere 0x10

8

128.2

f(Onset.Delta) 127.8

127.4

127.0 −1000 −1312 −183 −271 −376 −474 −573 −668 −792 −999

Onset.Delta

i.var="last.slope.weight", n.trees=bstop)

plot(fit, 9

180

160

f(last.slope.weight) 140

120

100

80 −5 0 5 10 15 20

last.slope.weight

i.var="alsfrs.score.slope", n.trees=bstop)

plot(fit,

## Warning in axis(side = 1, at = 1:227, labels = c("", "\025",

## "-0.334468864468864", : larghezza font sconosciuta per il carattere 0x15

10

180

160

f(alsfrs.score.slope) 140

120

100

80 −0.760916666666667 −2.62385057471264 1.18584415584416

alsfrs.score.slope

Dal dataset ALS otteniamo il seguente MSE sul test dato dall’utilizzo di 3 tecniche a noi note:

11

12

13

Boosting in maniera generale n

ˆ 1

[0] P

f arg min loss(y , c). b

1) Inizializzo la stima di = funzione di perdita = Fissiamo = 0 e un

c i

i=1

n

B >

numero di iterazioni grande con 0.

b r

2) Facciamo crescere di 1 e ricaviamo , il gradiente negativo della funzione di perdita, valutato al valore

i

−

b 1: ˆ

[b−1]

f f

(x) = (x )

i

δ

−

r loss[y , f i ..., n

= (x)]| per = 1,

ˆ

i i [b−1]

f f

(x)= (x )

δf (x) i

14

g(x)

3) Stimo/fitto il modello sui dati: [b]

→

, x ..., , x ĝ

(r ), (r ) (x)

n n

1 1

4) Aggiorno: ˆ ˆ

[b] [b−1] [b]

f f νĝ

(x) = (x) + (x)

< ν <

Dove 0 1 è la step-length factor.

b B B).

5) iteriamo i passi da 2 a 4 finchè non sarà uguale a (b = g(x) f

La funzione di perdita ci definisce il problema di minimo e il modello ci dice la forma funzionale che

può assumere.

Perdita di errore quadratico: 1 2

− f

loss[y, f [y (x)]

(x)] = 2

Problema di minimizzazione: 2

{[Y − } |X

arg min E f f E(Y x)

(X)] ha soluzione (x) = =

f Y,X

Abbiamo una funzione obiettivo che può essere qualsiasi funzione di perdita arbitraria.

Per semplicità abbiamo diviso la perdita quadratica per 2.

Il gradiente negativo è: δ 1 2

− − −

f y f

[y (x)] = (x)

δf (x) 2

dato dalla derivata della funzione di perdita con cambio di segno. ho diviso per 2 al fine di evitare di

raddoppiare il residuo.

Ci viene dato il boosting con un algoritmo di perdita di errore quadratico.

Perdita di Laplace (in valore assoluto): |y −

loss[y, f f

(x)] = (x)|

Problema di minimizzazione: {|Y − |X

arg min E f f M edian(Y x)

(X)|} ha soluzione (x) = =

f Y,X

Posso usare la funzione di perdita in valore assoluto (Laplace), ho un problema di minimo differente. La

soluzione in questo caso è la mediana che ci da risultati che vengono influenzati meno rispetto ad una

regressione lineare. Il gradiente negativo è:

δ

− |y − −

f sign(y f

(x)| = (x))

δf (x) y f

Da notare che la perdita di Laplace non è differenziabile in = (x); il valore del gradiente negativo in questo

punto è impostato sullo 0.

E’ l’approccio della regressione mediana. Abbiamo resistenza a coda lunga della distribuzione degli errori e

valori anomali (robustezza).

Sono libero di scegliere qualsiasi funzione di perdita. 15

Discesa del gradiente:

Immagino di costruire delle curve di livello. Questo algoritmo ci fa partire dal punto di partenza e ci fa

muovere in direzione ortogonale con un passo predefinito che ci porta nella funzione che ha il minimo, quindi

quando il gradiente scende di più. Devo stare attento a funzioni che hanno minimi globali, dato che rischio di

non arrivare al minimio globale, ma semplicemente a quello locale.

Se passo ad un problema di classificazione?

Classificazione ∈ {0,

Y

La variabile risposta è 1}. −1

Utilizzo spesso un’altra codifica in Machine Learning: 1 e dato che

− ∈ {−1,

Ỹ = 2Y 1 1}

− −1|X

π(x) P r(

Ỹ x), π(x) P r(

Ỹ x)

= = 1|X = 1 = = =

π(x)

1

f log

(x) = la metà della funzione logit

− π(x)

2 1

Ci riconduciamo al classificatore di Bayes:

 ⇔

se π(x) > f >

1 0.5 (x) 0





∗ −1 ⇔

C (x) = se π(x) < f <

0.5 (x) 0

 ⇔

se π(x) f

? = 0.5 (x) = 0



Ad ogni passo dell’algoritmo abbiamo la classificazione e la probabilità.

≤

ỹf

La missclassificazione avviene se e solo se (x) 0

Perdita esponenziale: loss[ỹ, f exp[−ỹf

(x)] = (x)]

Problema di minimizzazione:

π(x)

1

arg min E{exp[− Ỹ f f log

(x)]} ha soluzione (x) =

f − π(x)

2 1

Il gradiente negativo è: δ

− exp[−ỹf ỹ exp[−ỹf

(x)] = (x)]

δf (x) 16 ˆ

[b]

f

Quindi la stima boosting dopo b iterazioni corrisponde a (X) e classifica:

( ˆ

[b]

se f >

1 (x) 0

∗

C (x) = ˆ

[b]

−1 se f <

(x) 0

E le stime di probabilità sono: ˆ

[b]

exp[ f (x)]

[b]

π̂ (x) = ˆ ˆ

[b