Hands-On Algorithms for Predictive Modeling

Arthur Charpentier, Algorithms for Predictive Modeling - 2019
Hands-On Algorithms for Predictive Modeling*
(Penalized) Logistic Regression, Boosting, Neural Nets, SVM, Trees, Forests
Arthur Charpentier (Université du Québec à Montréal, Canada)
Third International Congress on Actuarial Science and Quantitative Finance
Manizales, Colombia, June 2019
@freakonometrics freakonometrics freakonometrics.hypotheses.org 1

Modeling a 0/1 random variable
Myocardial infarction of patients admited in E.R.
◦ heart rate (FRCAR) : x1
◦ heart index (INCAR) : x2
◦ stroke index (INSYS) : x3
◦ diastolic pressure (PRDIA) : x4
◦ pulmonary arterial pressure (PAPUL) : x5
◦ ventricular pressure (PVENT) : x6
◦ lung resistance (REPUL) : x7
◦ death or survival (PRONO) : y
dataset {(xi, yi)} where i = 1, . . . , n, n = 73.
1 > myocarde=read.table("http:// freakonometrics .free.fr/myocarde.csv",
head=TRUE ,sep=";")

Modeling a 0/1 random variable
1 > myocarde = read.table("http:// freakonometrics .free.fr/myocarde.csv"
,head=TRUE , sep=";")
2 > head(myocarde)
3 FRCAR INCAR INSYS PRDIA PAPUL PVENT REPUL PRONO
4 1 90 1.71 19.0 16 19.5 16.0 912 SURVIE
5 2 90 1.68 18.7 24 31.0 14.0 1476 DECES
6 3 120 1.40 11.7 23 29.0 8.0 1657 DECES
7 4 82 1.79 21.8 14 17.5 10.0 782 SURVIE
8 5 80 1.58 19.7 21 28.0 18.5 1418 DECES
9 6 80 1.13 14.1 18 23.5 9.0 1664 DECES
10 > myocarde$PRONO = (myocarde$PRONO =="SURVIE")*1

Classiﬁcation : Logistic Regression (and GLM)
Simulated data, y ∈ {0, 1}
P(Y = 1|X = x) =
exp[xTβ]
1 + exp[xTβ]
Inference using maximum likelihood techniques
β = argmin
n
i=1
log[P(Y = 1|X = x)]
and the score model is then
s(x) =
exp[xTβ]
1 + exp[xTβ]
1 x1 = c(.4 ,.55 ,.65 ,.9 ,.1 ,.35 ,.5 ,.15 ,.2 ,.85)
2 x2 = c(.85 ,.95 ,.8 ,.87 ,.5 ,.55 ,.5 ,.2 ,.1 ,.3)
3 y = c(1,1,1,1,1,0,0,1,0,0)
4 df = data.frame(x1=x1 ,x2=x2 ,y=as.factor(y))
q
q
q
q
q
q
q
q
q
q
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
q
q
q
q
q
q
q
q
q
q

Classiﬁcation : Logistic Regression (and GLM)
The logistic regression is based on the assumption that
Y |X = x ∼ B(px), where px =
exp[xTβ]
1 + exp[xTβ]
= 1 + exp[−xT
β]
−1
The goal is to estimate parameter β.
Recall that the heuristics for the use of that function for the probability is that
log[odds(Y = 1)] = log
P[Y = 1]
P[Y = 0]
= xT
β
• Maximum of the (log)-likelihood function
The log-likelihood is here
log L =
n
i=1
yi log pi + (1 − yi) log(1 − pi)
where pi = (1 + exp[−xi
Tβ])−1.

Mathematical Statistics Courses in a Nutshell
Consider observations {y1, · · · , yn}
from iid random variables Yi ∼ Fθ
(with “density” fθ).
Likelihood is
L(θ ; y) →
n
i=1
fθ(yi)
Maximum likelihood estimate is
θ
mle
∈ arg max
θ∈Θ
{L(θ; y)}
or θ
mle
∈ arg max
θ∈Θ
{log L(θ; y)}
Fisher (1912, On an absolute criterion for ﬁtting frequency curves).

Mathematical Statistics Courses in a Nutshell
Under standard assumptions (Identiﬁcation of the model, Compactness, Continuity and
Dominance), the maximum likelihood estimator is consistent θ
mle P
−→ θ. With additional
assumptions, it can be shown that the maximum likelihood estimator converges to a normal
distribution
√
n θ
mle
− θ
L
−→N(0, I−1
)
where I is Fisher information matrix (i.e. θ
mle
is asymptotically eﬃcient).
Eg. if Y ∼ N(θ, σ2), log L ∝ −
n
i=1
yi − θ
σ
2
, and θmle = y (see also method of moments).
max log L ⇐⇒ min
n
i=1
ε2
i = least squares

Classiﬁcation : Logistic Regression
• Numerical Optimization
Using optim() might not be a great idea.
BFGS algorithm (Broyden, Fletcher, Goldfarb & Shanno) based on Newton’s algorithm.
To ﬁnd a zero of function g, g(x0) = 0, use Taylor’s expansion,
0 = g(x0) = g(x) + g (x)(x0 − x)
i.e. (if x → xold and x0 → xnew)
xnew = xold −
g(xold)
g (xold)
from some starting point
or, in higher dimension
xnew = xold − g(xold)−1
· g(xold) from some starting point
Here g is the gradient of the log-likelihood,
βnew = βold −
∂2 log L(βold)
∂β∂βT
−1
·
∂ log L(βold)
∂β

Numerically, since only g is given, consider some small h > 0 and set
xnew = xold −
2h · g(xold)
g(xold + h) + g(xold − h)
from some starting point
1 y = myocarde$PRONO
2 X = cbind (1,as.matrix(myocarde [ ,1:7]))
3 negLogLik = function(beta){
4 -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta)))
5 }
6 beta_init = lm(PRONO˜.,data=myocarde)$ coefficients
7 logistic_opt = optim(par = beta_init , negLogLik , hessian=TRUE , method
= "BFGS", control=list(abstol =1e -9))

Here, we obtain
1 logistic_opt$par
2 (Intercept) FRCAR INCAR INSYS
3 1.656926397 0.045234029 -2.119441743 0.204023835
4 PRDIA PAPUL PVENT REPUL
5 -0.102420095 0.165823647 -0.081047525 -0.005992238
6 simu = function(i){
7 logistic_opt_i = optim(par = rnorm (8,0,3)*beta_init ,
8 negLogLik , hessian=TRUE , method = "BFGS",
9 control=list(abstol =1e-9))
10 logistic_opt_i$par [2:3]
11 }
12 v_beta = t(Vectorize(simu)(1:1000))

1 plot(v_beta)
2 par(mfrow=c(1 ,2))
3 hist(v_beta [,1], xlab=names(myocarde)[1])
4 hist(v_beta [,2], xlab=names(myocarde)[2])

• Fisher Algorithm
We want to use Newton’s algorithm
βnew = βold −
∂2 log L(βold)
∂β∂βT
−1
·
∂ log L(βold)
∂β
but here, the gradient and the hessian matrix have explicit expressions
∂ log L(βold)
∂β
= XT
(y − pold)
∂2 log L(βold)
∂β∂βT
= −XT
∆oldX
Note that it is fast and robust (to the starting point)
1 Y=myocarde$PRONO
2 X=cbind (1,as.matrix(myocarde [ ,1:7]))
3 colnames(X)=c("Inter",names(myocarde [ ,1:7]))
4 beta=as.matrix(lm(Y˜0+X)$coefficients ,ncol =1)
5 for(s in 1:9){
6 pi=exp(X%*%beta[,s])/(1+ exp(X%*%beta[,s]))
7 gradient=t(X)%*%(Y-pi)

8 omega=matrix (0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))
9 Hessian=-t(X)%*%omega%*%X
10 beta=cbind(beta ,beta[,s]-solve(Hessian)%*%gradient)}
11 beta [ ,8:10]
12 [,1] [,2] [,3]
13 XInter -10.187641685 -10.187641696 -10.187641696
14 XFRCAR 0.138178119 0.138178119 0.138178119
15 XINCAR -5.862429035 -5.862429037 -5.862429037
16 XINSYS 0.717084018 0.717084018 0.717084018
17 XPRDIA -0.073668171 -0.073668171 -0.073668171
18 XPAPUL 0.016756506 0.016756506 0.016756506
19 XPVENT -0.106776012 -0.106776012 -0.106776012
20 XREPUL -0.003154187 -0.003154187 -0.003154187

• iteratively reweighted least squares
We want to compute something like
βnew = (XT
∆oldX)−1
XT
∆oldz
(if we do substitute matrices in the analytical expressions) where
z = Xβold + ∆−1
old
[y − pold]
But actually, that’s simply a standard least-square problem
βnew = argmin (z − Xβ)T
∆−1
old
(z − Xβ)
The only problem here is that weights ∆old are functions of unknown βold
But actually, if we keep iterating, we should be able to solve it : given the β we got the weights,
and with the weights, we can use weighted OLS to get an updated β. That’s the idea of
iteratively reweighted least squares.
1 df = myocarde
2 beta_init = lm(PRONO˜.,data=df)$ coefficients
3 X = cbind (1,as.matrix(myocarde [ ,1:7]))
4 beta = beta_init
5 for(s in 1:1000){

6 p = exp(X %*% beta) / (1+ exp(X %*% beta))
7 omega = diag(nrow(df))
8 diag(omega) = (p*(1-p))
9 df$Z = X %*% beta + solve(omega) %*% (df$PRONO - p)
10 beta = lm(Z˜.,data=df[,-8], weights=diag(omega))$ coefficients
11 }
12 beta
13 (Intercept) FRCAR INCAR INSYS PRDIA
14 -10.187641696 0.138178119 -5.862429037 0.717084018 -0.073668171
15 PAPUL PVENT REPUL
16 0.016756506 -0.106776012 -0.003154187

1 summary( lm(Z˜.,data=df[,-8], weights=diag(omega)))
2
3 Coefficients :
4 Estimate Std. Error t value Pr(>|t|)
5 (Intercept) -10.187642 10.668138 -0.955 0.343
6 FRCAR 0.138178 0.102340 1.350 0.182
7 INCAR -5.862429 6.052560 -0.969 0.336
8 INSYS 0.717084 0.503527 1.424 0.159
9 PRDIA -0.073668 0.261549 -0.282 0.779
10 PAPUL 0.016757 0.306666 0.055 0.957
11 PVENT -0.106776 0.099145 -1.077 0.286
12 REPUL -0.003154 0.004386 -0.719 0.475

1 summary(glm(PRONO˜.,data=myocarde ,family=binomial(link = "logit")))
2
3 Coefficients :
4 Estimate Std. Error z value Pr(>|z|)
5 (Intercept) -10.187642 11.895227 -0.856 0.392
6 FRCAR 0.138178 0.114112 1.211 0.226
7 INCAR -5.862429 6.748785 -0.869 0.385
8 INSYS 0.717084 0.561445 1.277 0.202
9 PRDIA -0.073668 0.291636 -0.253 0.801
10 PAPUL 0.016757 0.341942 0.049 0.961
11 PVENT -0.106776 0.110550 -0.966 0.334
12 REPUL -0.003154 0.004891 -0.645 0.519

Classiﬁcation : Logistic Regression with Splines (and GAM)
To get home-made linear splines, use
β1x + β2(x − s)+
1 pos = function(x,s) (x-s)*(x>=s)
1 reg = glm(PRONO˜INSYS+pos(INSYS ,15)+pos(INSYS ,25) ,data=myocarde ,
family=binomial)
1 summary(reg)
2
3 Coefficients :
5 (Intercept) -0.1109 3.278 -0.034 0.973
6 INSYS -0.1751 0.252 -0.693 0.488
7 pos(INS ,15) 0.7900 0.374 2.109 0.034 *
8 pos(INS ,25) -0.5797 0.290 -1.997 0.045 *

Classiﬁcation : Logistic Regression with Splines
We can deﬁne spline functions with support (5, 55) and with knots {15, 25}
1 library(splines)
2 clr6 = c("#1 b9e77","#d95f02","#7570 b3","#
e7298a","#66 a61e","#e6ab02")
3 x = seq (0,60,by =.25)
4 B = bs(x,knots=c(15 ,25) ,Boundary.knots=c(5 ,55)
,degre =1)
5 matplot(x,B,type="l",lty=1,lwd=2,col=clr6)

1 reg = glm(PRONO˜bs(INSYS ,knots=c(15 ,25) ,
2 Boundary.knots=c(5 ,55) ,degre =1) ,
3 data=myocarde ,family=binomial)
4 summary(reg)
5
6 Coefficients :
8 (Intercept) -0.9863 2.055 -0.480 0.631
9 bs(INSYS)1 -1.7507 2.526 -0.693 0.488
10 bs(INSYS)2 4.3989 2.061 2.133 0.032 *
11 bs(INSYS)3 5.4572 5.414 1.008 0.313

For quadratic splines, consider a decomposition on x,x2,(x − s1)2
+
1 pos2 = function(x,s) (x-s)ˆ2*(x>=s)
2 reg = glm(PRONO˜poly(INSYS ,2)+pos2(INSYS ,15)+
pos2(INSYS ,25) ,
3 data=myocarde ,family=binomial)
4 summary(reg)
5
6 Coefficients :
8 (Intercept) 29.9842 15.236 1.968 0.049 *
9 poly (2)1 408.7851 202.419 2.019 0.043 *
10 poly (2)2 199.1628 101.589 1.960 0.049 *
11 pos2 (15) -0.2281 0.126 -1.805 0.071 .
12 pos2 (25) 0.0439 0.080 0.545 0.585

1 x = seq (0,60,by =.25)
2 B=bs(x,knots=c(15 ,25) ,Boundary.knots=c(5 ,55) ,
degre =2)
3 matplot(x,B,type="l",xlab="INSYS",col=clr6)

1 reg = glm(PRONO˜bs(INSYS ,knots=c(15 ,25) ,
2 Boundary.knots=c(5 ,55) ,degre =2) ,data=myocarde ,
3 family=binomial)
4 summary(reg)
5
6 Coefficients :
8 (Intercept) 7.186 5.261 1.366 0.172
9 bs(INSYS)1 -14.656 7.923 -1.850 0.064 .
10 bs(INSYS)2 -5.692 4.638 -1.227 0.219
11 bs(INSYS)3 -2.454 8.780 -0.279 0.779
12 bs(INSYS)4 6.429 41.675 0.154 0.877
Finally, for cubic splines, use x, x2, x3, (x − s1)3
+, (x − s2)3
+, etc.

1 B=bs(x,knots=c(15 ,25) ,Boundary.knots=c(5 ,55) ,
degre =3)
2 matplot(x,B,type="l",lwd=2,col=clr6 ,lty=1,ylim
=c( -.2 ,1.2))
3 abline(v=c(5 ,15 ,25 ,55) ,lty =2)
1 reg = glm(PRONO˜1+bs(INSYS ,degree =1,df =4),data
=myocarde ,family=binomial)
2 attr(reg$terms , "predvars")[[3]]
3 bs(INSYS , degree = 1L, knots = c(15.8 , 21.4 ,
27.15) ,
4 Boundary.knots = c(8.7 , 54) , intercept = FALSE
)
5 quantile(myocarde$INSYS ,(0:4)/4)
6 0% 25% 50% 75% 100%
7 8.70 15.80 21.40 27.15 54.00

1 B = bs(x,degree =1,df=4)
2 B = cbind (1,B)
3 y = B%*% coefficients (reg)
4 plot(x,y,type="l",col="red",lwd =2)
5 abline(v=quantile(myocarde$INSYS ,(0:4)/4),lty
=2)

1 reg = glm(y˜bs(x1 ,degree =1,df =3)+bs(x2 ,degree
=1,df =3) ,data=df ,family=binomial(link = "
logit"))
2 u = seq(0,1, length =101)
3 p = function(x,y) predict.glm(reg ,newdata=data
.frame(x1=x,x2=y),type="response")
4 v = outer(u,u,p)
5 image(u,u,v,xlab="Variable 1",ylab="Variable 2
",col=clr10 ,breaks =(0:10)/10)
6 points(df$x1 ,df$x2 ,pch =19, cex =1.5 , col="white")
7 points(df$x1 ,df$x2 ,pch=c(1 ,19) [1+( df$y=="1")],
cex =1.5)
8 contour(u,u,v,levels = .5,add=TRUE)

Classiﬁcation : Regression with Kernels and k-NN
We do not use the probit framework, it is purely a non-parametric one.
The moving regressogram is
ˆm(x) =
n
i=1
yi · 1(xi ∈ [x ± h])
n
i=1
1(xi ∈ [x ± h])
Consider
˜m(x) =
n
i=1
yi · kh (x − xi)
n
i=1
kh (x − xi)
where (classically), for some kernel k,
kh(x) = k
x
h
k
−2 −1 0 1 2

1 mean_x = function(x,bw){
2 w = dnorm (( myocarde$INSYS -x)/bw , mean=0,sd
=1)
3 weighted.mean(myocarde$PRONO ,w)}
4 u = seq (5,55, length =201)
5 v = Vectorize(function(x) mean_x(x,3))(u)
6 plot(u,v,ylim =0:1 , type="l",col="red")
7 points(myocarde$INSYS ,myocarde$PRONO ,pch =19)
1 v = Vectorize(function(x) mean_x(x,2))(u)
2 plot(u,v,ylim =0:1 , type="l",col="red")

One can also use the ksmooth() function of R,
1 reg = ksmooth(myocarde$INSYS ,myocarde$PRONO ,"
normal",bandwidth = 2*exp (1))
2 plot(reg$x,reg$y,ylim =0:1 , type="l",col="red",
lwd=2,xlab="INSYS",ylab="")
Classical bias / variance tradeoﬀ (below mh(15) and mh(25))
10 20 30 40 50
0.00.20.40.60.81.0
INSYS
q
qq
q
qq
qq qq qq q qqqqq qq q
q qq qqq q
qqq
q q
q
q
qqq
q
q
qq
qq
q
q
q q
q
q
q q
q q
q
q
q q q q
q
q
q
q qq
qq
q
q
q
Density
0.0 0.2 0.4 0.6 0.8 1.0
0123456
Density
0.0 0.2 0.4 0.6 0.8 1.0
0123456

1 u = seq(0,1, length =101)
2 p = function(x,y){
3 bw1 = .2; bw2 = .2
4 w = dnorm ((df$x1 -x)/bw1 , mean=0,sd=1)*
5 dnorm ((df$x2 -y)/bw2 , mean=0,sd =1)
6 weighted.mean(df$y=="1",w)
7 }
8 v = outer(u,u,Vectorize(p))
9 image(u,u,v,col=clr10 ,breaks =(0:10)/10)
cex =1.5)

One can also use k-nearest neigthbors
˜mk(x) =
1
n
n
i=1
ωi,k(x)yi
where ωi,k(x) = n/k if i ∈ Ik
x with
Ik
x = {i : xi one of the k nearest observations to x}I
Classically, use Mahalanobis distance
1 Sigma = var(myocarde [ ,1:7])
2 Sigma_Inv = solve(Sigma)
3 d2_ mahalanobis = function(x,y,Sinv){as.numeric(x-y)%*%Sinv%*%t(x-y)}
4 k_closest = function(i,k){
5 vect_dist = function(j) d2_ mahalanobis (myocarde[i ,1:7] , myocarde[j
,1:7] , Sigma_Inv)
6 vect = Vectorize(vect_dist)((1: nrow(myocarde)))
7 which (( rank(vect)))}

1 k_majority = function(k){
2 Y=rep(NA ,nrow(myocarde))
3 for(i in 1: length(Y)) Y[i] = sort(myocarde$PRONO[k_closest(i,k)])[(
k+1)/2]
4 return(Y)}
1 k_mean = function(k){
2 Y=rep(NA ,nrow(myocarde))
3 for(i in 1: length(Y)) Y[i] = mean(myocarde$PRONO[k_closest(i,k)])
4 return(Y)}
1 Sigma_Inv = solve(var(df[,c("x1","x2")]))
2 u = seq(0,1, length =51)
4 k = 6
5 vect_dist = function(j) d2_ mahalanobis(c(x,y),df[j,c("x1","x2")],
Sigma_Inv)
6 vect = Vectorize(vect_dist)(1: nrow(df))
7 idx = which(rank(vect) <=k)

8 return(mean ((df$y==1)[idx ]))}
1 v = outer(u,u,Vectorize(p))
2 image(u,u,v,col=clr10 ,breaks =(0:10)/10)
cex =1.5)

Goodness of Fit: ROC Curve
Confusion matrix
Given a sample (yi, xi) and a model m, the confusion matrix is the contin-
gency table, with dimensions observed yi ∈ {0, 1} and predicted yi ∈ {0, 1}.
y
0(-) 1(+)
0(-) TN FN
y
1(+) FP TP
FP
TN+FP
TP
FN+TP
FPR TPR
Classical measures are
true positive rate (TPR) - or sensitivity
false positive rate (FPR) - or fall-out,
true negative rate (TNR) - or speciﬁcity
TNR = 1-FPR
among others (see wikipedia)
ROCR::performance(prediction(Score,Y),"tpr","fpr")

ROC (Receiver Operating Characteristic ) Curve
Assume that mt is deﬁne from a score function s, with mt(x) = 1(s(x) > t)
for some threshold t. The ROC curve is the curve (FPRt, TPRt obtained
from confusion matrices of mt’s.
n = 100 individuals
50 yi = 0 and 50 yi = 1
(well balanced)
q
q
Y=0 Y=1
OBSERVED
Y=1Y=0
PREDICTED
25 25
25 25
FPR TPR
25 25
50 50
~ 0.5 ~ 0.5
q
q
FALSE POSITIVE RATE
TRUEPOSITIVERATE
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
q

29 deaths (y = 0) and 42 survivals (y = 1)
y =
1 if P[Y = 1|X] > 0%
0 if P[Y = 1|X] ≤ 0%
0%
y
0 1
0 0 0
y
1 29 42
29
29+0
42
42+0
= 100% = 100%
FPR TPR

y =
1 if P[Y = 1|X] > 15%
0 if P[Y = 1|X] ≤ 15%
15%
y
0 1
0 17 2
y
1 12 40
12
17+12
40
42+2
∼ 41.4% ∼ 95.2%
FPR TPR

y =
1 if P[Y = 1|X] > 50%
0 if P[Y = 1|X] ≤ 50%
50%
y
0 1
0 25 3
y
1 4 39
4
25+4
39
39+3
∼ 13.8% ∼ 92.8%
FPR TPR

y =
1 if P[Y = 1|X] > 85%
0 if P[Y = 1|X] ≤ 85%
50%
y
0 1
0 28 13
y
1 1 29
1
28+1
29
29+13
∼ 3.4% ∼ 69.9%
FPR TPR

y =
1 if P[Y = 1|X] > 100%
0 if P[Y = 1|X] ≤ 100%
50%
y
0 1
0 29 42
y
0 0 29
0
29+0
0
42+0
= 0.0% = 0.0%
FPR TPR

Goodness of Fit
Be carefull of overfit...
Important to use
training / validation
• classification tree
• boosting classifier

ROC (Receiver Operating Characteristic ) Curve
Assume that mt is deﬁne from a score function s, with mt(x) = 1(s(x) > t)
for some threshold t. The ROC curve is the curve (FPRt, TPRt obtained
from confusion matrices of mt’s.
qq qqq qq
q
qq qqqq qqqq
q
qqqqq
q
qqq qqq qq
qq
qqqqqqqqqq qqq qq
q
qq qqq qqqq qq qq qq qqq
q
q qq qqq q qq qq q qqqqq qqq qqqq qqqq qq qqqq q
q
qqq qq qqq qqq qqq qq qqq qqqqqqq q qqqqq q qqq
q
qq qq qqq
q
q
q
qqq qq qqqq q qq q qqqq qqq q qq qqq qqq qq qqqq qqqq
q
q qqq qqqq
q
qq
q
q
qq
qq qqq q qqqq
q
qqq qqq qqqqqq qq q
q
qq qqq qqq q qqqqq qqq q
q
q qq qq qqqq
q
qqqq qqqqq q qqq qq
q
qqq q
q
q q qqq
q
qq qqq qqq qq qqq q
q
q qq qqq qq
q
q qqqq qq
q
q q
q
qq q
qq
qqq
q
q qqq q
q
qqqq qqqqq qqqq qqqq qqq
q
q qq qq qqqq q
q
qqq
q
qqq qqqq qq q qq q qq qqqqqq q qqqq q qqqq qqqqq
q
qq q qq qqqq qqqq
q
q
qq
qq q
q
qqq qqqq q
q
q
q
qq qq
qq
qqq q qq q qqqq qqqq qqqq qq qq qq qqqq qqqq q qqqq qqqq
q
q qqqq qqq qqq q
q
qq qqqq qqqqqqqq qq qqqq qqq q
q
qq qq
q
q q qqqq qqq q qqq q q qqqqqq q
q
q q qqqqq q q qqqq qqqq q qqqq qqqqq qq qqq q
q
qqq qq qqqq qq qq
qq
qqqqqq qq qqq
q
q qqqqqq
q
q qq qq q qqq q qqq q qq qq qqq q
q
qq q
q
q qq qq qq
q
q qqq qqq qqqqq qqqq q qqqqq
q
q
q
q q
qq
qq
q
qqqq
q
q qq q q
q
q
q
qqqq qq qq q qqq q qq qqqqqq qq qq qqqq qqqq q qqqq qq q qq qq qqqq q qq q
q
qqqqq qqqq qq qq qqqq q qqq
q
qq
q
qq qqqq qqqqq qqqq qq qqq qqq qqqq qq q q
q
qqq qq qq
q
qq qq qqq qq qqqqq qqq qqq qqqqq qqqqq qq qqq qqqq qqq q qq qq qqqq q qq qq qqqqq qq qqqq
qq
qqqqqqqq qqqq qq q
q
qqqqq qqq qq
q
qqqqq
q
q q
q
qqq qq q qqqq qq q qqqq qqqq qqqqq qq qqqqq
q
q q qq qqq qqq q
q
q
q
q qq q qqqq
q
qqq q qq qq q
SCORE S
OBSERVATIONY
0.0 0.2 0.4 0.6 0.8 1.0
Y=0
Y=1
FALSE POSITIVE RATE
TRUEPOSITIVERATE
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
q

With categorical variables, we have a collection of points
Need to interpolate between those points to have a curve, see convexiﬁcation of ROC curves,
e.g. Flach (2012, Machine Learning)
q
q
FALSE SPAM RATE
TRUESPAMRATE
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
q
q
q
q
q
q
FALSE SPAM RATE
TRUESPAMRATE
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
q
q

Goodness of Fit
One can derive a conﬁdence interval of the ROC curve using
boostrap techniques,
pROC::ci.se()
Various measures can be used see library hmeasures , with Gini
index, the Area Under the curve, or Kolmogorov-Smirnov
ks = sup
t∈R
|F1(t) − F0(t)|
ks = sup
t∈R
|
1
n1
i:yi=1
1s(xi)≤t −
1
n0
i:yi=0
1s(xi)≤t|
Specificity (%)
Sensitivity(%)
020406080100
100 80 60 40 20 0
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
Fn(x)
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
Survival
Death

Goodness of Fit
The Area Under the Curve, AUC, can be interpreted as the probability that a classiﬁer will
rank a randomly chosen positive instance higher than a randomly chosen negative one, see
Swets, Dawes & Monahan (2000, Psychological Science Can Improve Diagnostic Decisions)
Kappa statistic κ compares an Observed Accuracy with an Expected Accuracy (random
chance), see Landis & Koch (1977, The Measurement of Observer Agreement for Categorical Data
)
Y = 0 Y = 1
Y = 0 TN FN TN+FN
Y = 1 FP TP FP+TP
TN+FP FN+TP n
See also Obsersed and Random Confusion Tables
Y = 0 Y = 1
Y = 0 25 3 28
Y = 1 4 39 43
29 42 71
Y = 0 Y = 1
Y = 0 11.44 16.56 28
Y = 1 17.56 25.44 43
29 42 71

Goodness of Fit
Y = 0 Y = 1
Y = 0 25 3 28
Y = 1 4 39 43
29 42 71
Y = 0 Y = 1
Y = 0 11.44 16.56 28
Y = 1 17.56 25.44 43
29 42 71
total accuracy =
TP + TN
n
∼ 90.14%
random accuracy =
[TN + FP] · [TP + FN] + [TP + FP] · [TN + FN]
n2
∼ 51.93%
κ =
total accuracy − random accuracy
1 − random accuracy
∼ 79.48%

Penalized Estimation
Consider a parametric model, with true (unkrown) parameter θ, then
mse(ˆθ) = E (ˆθ − θ)2
= E (ˆθ − E ˆθ )2
variance
+ E (E ˆθ − θ)2
bias2
One can think of a shrinkage of an unbiased estimator,
Let θ denote an unbiased estimator of θ. Then
ˆθ =
θ2
θ2 + mse(θ)
· θ = θ −
mse(θ)
θ2 + mse(θ)
· θ
penalty
satisﬁes mse(ˆθ) ≤ mse(θ).
−2 −1 0 1 2 3 4
0.00.20.40.60.8
variance

Linear Regression Shortcoming
Least Squares Estimator β = (XT
X)−1XT
y
Unbiased Estimator E[β] = β
Variance Var[β] = σ2(XT
X)−1
which can be (extremely) large when det[(XT
X)] ∼ 0.
X =





1 −1 2
1 0 1
1 2 −1
1 1 0





then XT
X =



4 2 2
2 6 −4
2 −4 6


 while XT
X + I =



5 2 2
2 7 −4
2 −4 7



eigenvalues : {10, 6, 0} {11, 7, 1}
Ad-hoc strategy: use XT
X + λI

Linear Regression Shortcoming
Evolution of (β1, β2) →
n
i=1
[yi − (β1x1,i + β2x2,i)]2
when cor(X1, X2) = r ∈ [0, 1], on top.
Below, Ridge regression
(β1, β2) →
n
i=1
[yi − (β1x1,i + β2x2,i)]2
+λ(β2
1 + β2
2 )
where λ ∈ [0, ∞), below,
when cor(X1, X2) ∼ 1 (colinearity).
−2 −1 0 1 2 3 4
−3−2−10123
β1
β2
500
1000
1500
2000
2000
2500
2500 2500
2500
3000
3000 3000
3000
3500
q
−2 −1 0 1 2 3 4
−3−2−10123
β1
β2
1000
1000
2000
2000
3000
3000
4000
4000
5000
5000
6000
6000
7000

Normalization : Euclidean 2 vs. Mahalonobis
We want to penalize complicated models :
if βk is “too small”, we prefer to have βk = 0.
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Instead of d(x, y) = (x − y)T
(x − y)
use dΣ(x, y) = (x − y)TΣ−1(x − y)
beta1
beta2
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

Ridge Regression
... like the least square, but it shrinks estimated coeﬃcients towards 0.
β
ridge
λ = argmin
n
i=1
(yi − xT
i β)2
+ λ
p
j=1
β2
j
β
ridge
λ = argmin



y − Xβ
2
2
=criteria
+ λ β 2
2
=penalty



λ ≥ 0 is a tuning parameter.
The constant is usually unpenalized. The true equation is
β
ridge
λ = argmin



y − (β0 + Xβ)
2
2
=criteria
+ λ β
2
2
=penalty




Ridge Regression
β
ridge
λ = argmin y − (β0 + Xβ)
2
2
+ λ β
2
2
can be seen as a constrained optimization problem
β
ridge
λ = argmin
β 2
2
≤hλ
y − (β0 + Xβ)
2
2
Explicit solution
βλ = (XT
X + λI)−1
XT
y
If λ → 0, β
ridge
0 = β
ols
If λ → ∞, β
ridge
∞ = 0.
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
30
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

Ridge Regression
This penalty can be seen as rather unfair if components of
x are not expressed on the same scale
• center: xj = 0, then β0 = y
• scale: xT
j xj = 1
Then compute
β
ridge
λ = argmin



y − Xβ 2
2
=loss
+ λ β 2
2
=penalty



beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
30
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
40
40
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

Ridge Regression
Observe that if xj1 ⊥ xj2 , then
β
ridge
λ = [1 + λ]−1
β
ols
λ
which explain relationship with shrinkage.
But generally, it is not the case...
q
q
Theorem There exists λ such that mse[β
ridge
λ ] ≤ mse[β
ols
λ ]

Ridge Regression
Lλ(β) =
n
i=1
(yi − β0 − xT
i β)2
+ λ
p
j=1
β2
j
∂Lλ(β)
∂β
= −2XT
y + 2(XT
X + λI)β
∂2Lλ(β)
∂β∂βT
= 2(XT
X + λI)
where XT
X is a semi-positive deﬁnite matrix, and λI is a positive deﬁnite matrix, and
βλ = (XT
X + λI)−1
XT
y

The Bayesian Interpretation
From a Bayesian perspective,
P[θ|y]
posterior
∝ P[y|θ]
likelihood
· P[θ]
prior
i.e. log P[θ|y] = log P[y|θ]
log likelihood
+ log P[θ]
penalty
If β has a prior N(0, τ2I) distribution, then its posterior distribution has mean
E[β|y, X] = XT
X +
σ2
τ2
I
−1
XT
y.

Properties of the Ridge Estimator
βλ = (XT
X + λI)−1
XT
y
E[βλ] = XT
X(λI + XT
X)−1
β.
i.e. E[βλ] = β.
Observe that E[βλ] → 0 as λ → ∞.
Assume that X is an orthogonal design matrix, i.e. XT
X = I, then
βλ = (1 + λ)−1
β
ols
.

Set W λ = (I + λ[XT
X]−1)−1. One can prove that
W λβ
ols
= βλ.
Thus,
Var[βλ] = W λVar[β
ols
]W T
λ
and
Var[βλ] = σ2
(XT
X + λI)−1
XT
X[(XT
X + λI)−1
]T
.
Observe that
Var[β
ols
] − Var[βλ] = σ2
W λ[2λ(XT
X)−2
+ λ2
(XT
X)−3
]W T
λ ≥ 0.

Hence, the conﬁdence ellipsoid of ridge estimator is indeed
smaller than the OLS,
If X is an orthogonal design matrix,
Var[βλ] = σ2
(1 + λ)−2
I.
mse[βλ] = σ2
trace(W λ(XT
X)−1
W T
λ) + βT
(W λ − I)T
(W λ − I)β.
If X is an orthogonal design matrix,
mse[βλ] =
pσ2
(1 + λ)2
+
λ2
(1 + λ)2
βT
β
0.0 0.2 0.4 0.6 0.8
−1.0−0.8−0.6−0.4−0.2
β1
β2
1
2
3
4
5
6
7

mse[βλ] =
pσ2
(1 + λ)2
+
λ2
(1 + λ)2
βT
β
is minimal for
λ =
pσ2
βT
β
Note that there exists λ > 0 such that mse[βλ] < mse[β0] = mse[β
ols
].

Sparsity Issues
In severall applications, k can be (very) large, but a lot of features are just noise: βj = 0 for
many j’s. Let s denote the number of relevent features, with s << k, cf Hastie, Tibshirani &
Wainwright (2015) Statistical Learning with Sparsity,
s = card{S} where S = {j; βj = 0}
The model is now y = XT
S βS + ε, where XT
S XS is a full rank matrix.
q
q
= . +

Going further on sparcity issues
The Ridge regression problem was to solve
β = argmin
β∈{ β 2
≤s}
{ Y − XT
β 2
2
}
Deﬁne a 0 = 1(|ai| > 0).
Here dim(β) = k but β 0 = s.
We wish we could solve
β = argmin
β∈{ β 0
=s}
{ Y − XT
β 2
2
}
Problem: it is usually not possible to describe all possible constraints, since
s
k
coeﬃcients
should be chosen here (with k (very) large).
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

In a convex problem, solve the dual problem,
e.g. in the Ridge regression : primal problem
min
β∈{ β 2
≤s}
{ Y − XT
β 2
2
}
and the dual problem
min
β∈{ Y −XTβ 2
≤t}
{ β 2
2
}
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
26
27
30
32
35
40
40
50
60
70
80
90
100
110
120
120
130
130
140 140
X
q
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
26
27
30
32
35
40
40
50
60
70
80
90
100
110
120
120
130
130
140 140
X
q
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

Idea: solve the dual problem
β = argmin
β∈{ Y −XTβ 2
≤h}
{ β 0 }
where we might convexify the 0 norm, · 0 .

On [−1, +1]k, the convex hull of β 0 is β 1
On [−a, +a]k, the convex hull of β 0 is a−1 β 1
Hence, why not solve
β = argmin
β; β 1
≤˜s
{ Y − XT
β 2 }
which is equivalent (Kuhn-Tucker theorem) to the Lagragian optimization problem
β = argmin{ Y − XT
β 2
2
+λ β 1 }

lasso Least Absolute Shrinkage and Selection Operator
β ∈ argmin{ Y − XT
β 2
2
+λ β 1 }
is a convex problem (several algorithms ), but not strictly convex (no unicity of the minimum).
Nevertheless, predictions y = xTβ are unique.
MM, minimize majorization, coordinate descent Hunter & Lange (2003, A Tutorial on MM
Algorithms).

lasso Regression
No explicit solution...
If λ → 0, β
lasso
0 = β
ols
If λ → ∞, β
lasso
∞ = 0.
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

lasso Regression
For some λ, there are k’s such that β
lasso
k,λ = 0.
Further, λ → β
lasso
k,λ is piecewise linear
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
30
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0
beta1
beta2
−1 −0.5 0.5 1
−1
−0.5
0.5
1
30
40
40
50
60
70
80
90
100
110
120
120
150 150
40
40
X
−1.0 −0.5 0.0 0.5 1.0
−1.0−0.50.00.51.0

lasso Regression
In the orthogonal case, XT
X = I,
β
lasso
k,λ = sign(β
ols
k ) |β
ols
k | −
λ
2
i.e. the LASSO estimate is related to the soft threshold
function...
q
q

Classiﬁcation : Penalized Logistic Regression, 2 norm (Ridge)
We want to solve something like
βλ = argmin{− log L(β|x, y) + λ β 2
2
}
In the case where we consider the log-likelihood of some Gaussian variable, we get the sum of
the square of the residuals, and we can obtain an explicit solution. But not in the context of a
logistic regression.
The heuristics about Ridge regression is the following graph. In the background, we can
visualize the (two-dimensional) log-likelihood of the logistic regression, and the blue circle is the
constraint we have, if we rewrite the optimization problem as a constrained optimization
problem :
min
β: β 2
2
≤s
n
i=1
− log L(yi, β0 + xT
β)

It can be written equivalently (it is a strictly convex problem)
min
β,λ
−
n
i=1
log L(yi, β0 + xT
β) + λ β 2
2
Thus, the constrained maximum should lie in the blue disk.
1 > PennegLogLik = function(bbeta ,lambda =0){
2 b0 = bbeta [1]
3 beta = bbeta [-1]
4 -sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*
5 log (1 + exp(b0+X%*%beta)))+lambda*sum(beta ˆ2)}

1 > opt_ridge = function(lambda){
2 beta_init = lm(PRONO˜.,data=myocarde)$
coefficients
3 logistic_opt = optim(par = beta_init*0,
function(x) PennegLogLik (x,lambda), method
= "BFGS", control=list(abstol =1e -9))
4 logistic_opt$par [ -1]}
5 > v_lambda = c(exp(seq(-2,5, length =61)))
6 > est_ridge = Vectorize(opt_ridge)(v_lambda)
7 > plot(v_lambda ,est_ridge [1 ,])

• Using Fisher’s (penalized) score
On the penalized problem, we can easily prove that
∂ log Lp(βλ,old)
∂β
=
∂ log L(βλ,old)
∂β
− 2λβold
while
∂2 log Lp(βλ,old)
∂β∂βT
=
∂2 log L(βλ,old)
∂β∂βT
− 2λI
Hence
βλ,new = (XT
∆oldX + 2λI)−1
XT
∆oldz
And interestingly, with that algorithm, we can also derive the variance of the estimator
Var[βλ] = [XT
∆X + 2λI]−1
XT
∆Var[z]∆X[XT
∆X + 2λI]−1
where Var[z] = ∆−1

The code to compute βλ is then
1 Y = myocarde$PRONO
2 X = myocarde [ ,1:7]
3 for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
4 X = as.matrix(X)
5 X = cbind (1,X)
6 colnames(X) = c("Inter",names(myocarde [ ,1:7]))
7 beta = as.matrix(lm(Y˜0+X)$coefficients ,ncol =1)
8 for(s in 1:9){
9 pi = exp(X%*%beta[,s])/(1+ exp(X%*%beta[,s]))
10 Delta = matrix (0,nrow(X),nrow(X));diag(Delta)=(pi*(1-pi))
11 z = X%*%beta[,s] + solve(Delta)%*%(Y-pi)
12 B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%
Delta%*%z)
13 beta = cbind(beta ,B)}
14 beta [ ,8:10]
15 [,1] [,2] [,3]

16 XInter 0.59619654 0.59619654 0.59619654
17 XFRCAR 0.09217848 0.09217848 0.09217848
18 XINCAR 0.77165707 0.77165707 0.77165707
19 XINSYS 0.69678521 0.69678521 0.69678521
20 XPRDIA -0.29575642 -0.29575642 -0.29575642
21 XPAPUL -0.23921101 -0.23921101 -0.23921101
22 XPVENT -0.33120792 -0.33120792 -0.33120792
23 XREPUL -0.84308972 -0.84308972 -0.84308972
Use
Var[βλ] = [XT
∆X + 2λI]−1
XT
∆Var[z]∆X[XT
∆X + 2λI]−1
where Var[z] = ∆−1.
1 newton_ridge = function(lambda =1){
2 beta = as.matrix(lm(Y˜0+X)$coefficients ,ncol =1)*runif (8)
3 for(s in 1:20){
4 pi = exp(X%*%beta[,s])/(1+ exp(X%*%beta[,s]))
5 Delta = matrix (0,nrow(X),nrow(X));diag(Delta)=(pi*(1-pi))
6 z = X%*%beta[,s] + solve(Delta)%*%(Y-pi)

7 B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%
Delta%*%z)
8 beta = cbind(beta ,B)}
9 Varz = solve(Delta)
10 Varb = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% t(X)%*%
Delta %*% Varz %*%
11 Delta %*% X %*% solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X)))
12 return(list(beta=beta[,ncol(beta)],sd=sqrt(diag(Varb))))}
1 v_lambda=c(exp(seq(-2,5, length =61)))
2 est_ridge=Vectorize(function(x) newton_ridge(x
)$beta)(v_lambda)
3 library(" RColorBrewer ")
4 colrs=brewer.pal(7,"Set1")
5 plot(v_lambda ,est_ridge [1,],col=colrs [1], type=
"l")
6 for(i in 2:7) lines(v_lambda ,est_ridge[i,],col
=colrs[i])

And for the variance
2 est_ridge=Vectorize(function(x) newton_ridge(x
)$sd)(v_lambda)
3 library(" RColorBrewer ")
4 colrs=brewer.pal(7,"Set1")
5 plot(v_lambda ,est_ridge [1,],col=colrs [1], type=
"l")
6 for(i in 2:7) lines(v_lambda ,est_ridge[i,],col
=colrs[i],lwd =2)

3 for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X
[,j])
4 X = as.matrix(X)
5 library(glmnet)
6 glm_ridge = glmnet(X, y, alpha =0)
7 plot(glm_ridge ,xvar="lambda",col=colrs ,lwd =2)

• Ridge Regression with Orthogonal Covariates
1 library(factoextra)
2 pca = princomp(X)
3 pca_X = get_pca_ind(pca)$coord
4
5 library(glmnet)
6 glm_ridge = glmnet(pca_X, y, alpha =0)
7 plot(glm_ridge ,xvar="lambda",col=colrs ,lwd =2)
8
9 plot(glm_ridge ,col=colrs ,lwd =2)

1 df0 = df
2 df0$y=as.numeric(df$y) -1
3 plot_lambda = function(lambda){
4 m = apply(df0 ,2,mean)
5 s = apply(df0 ,2,sd)
6 for(j in 1:2) df0[,j] = (df0[,j]-m[
j])/s[j]
7 reg = glmnet(cbind(df0$x1 ,df0$x2),
df0$y==1, alpha =0)
9 plot(reg ,xvar="lambda",col=c("blue"
,"red"),lwd =2)
10 abline(v=log (.2))
11 plot_lambda (.2)
12 plot(reg ,xvar="lambda",col=c("blue"
,"red"),lwd =2)
13 abline(v=log (1.2))
14 plot_lambda (1.2)

Classiﬁcation : Penalized Logistic Regression, 1 norm (lasso)
It can be written equivalently (it is a strictly convex problem)
min
β,λ
−
n
i=1
log L(yi, β0 + xT
β) + λ β 1
Thus, the constrained maximum should lie in the blue disk.

3 for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X
[,j])
4 X = as.matrix(X)
5 LogLik = function(bbeta){
6 b0=bbeta [1]
7 beta=bbeta [-1]
8 sum(-y*log(1 + exp(-(b0+X%*%beta))) -
9 (1-y)*log(1 + exp(b0+X%*%beta)))}
10 u = seq(-4,4, length =251)
11 v = outer(u,u,function(x,y) LogLik(c(1,x,y)))
12 image(u,u,v,col=rev(heat.colors (25)))
13 contour(u,u,v,add=TRUE)
14 polygon(c(-1,0,1,0),c(0,1,0,-1),border="blue")

1 PennegLogLik = function(bbeta ,lambda =0){
2 b0=bbeta [1]
3 beta=bbeta [-1]
4 -sum(-y*log(1 + exp(-(b0+X%*%beta))) -
5 (1-y)*log(1 + exp(b0+X%*%beta)))+lambda*sum(abs(beta)) }
6 opt_lasso = function(lambda){
7 beta_init = lm(PRONO˜.,data=myocarde)$ coefficients
8 logistic_opt = optim(par = beta_init*0, function(x) PennegLogLik (x,
lambda),
9 hessian=TRUE , method = "BFGS", control=list(abstol =1e -9))
1 logistic_opt$par [-1] }
3 est_lasso=Vectorize(opt_lasso)(v_lambda)
4 plot(v_lambda ,est_lasso [1,],col=colrs [1])
5 for(i in 2:7) lines(v_lambda ,est_lasso[i,],col
=colrs[i],lwd =2)

Observe that coeﬃcients are not (strictly) null
1 library(glmnet)
2 glm_lasso = glmnet(X, y, alpha =1)
3 plot(glm_lasso ,xvar="lambda",col=colrs ,lwd =2)
4 plot(glm_lasso ,col=colrs ,lwd =2)
5 glmnet(X, y, alpha =1, lambda=exp (-4))$beta
6 7x1 sparse Matrix of class "dgCMatrix"
7 s0
8 FRCAR .
9 INCAR 0.11005070
10 INSYS 0.03231929
11 PRDIA .
12 PAPUL .
13 PVENT -0.03138089
14 REPUL -0.20962611

With orthogonal covariates
4 library(factoextra)
5 pca = princomp(X)
6 pca_X = get_pca_ind(pca)$coord
7 glm_lasso = glmnet(pca_X, y,
alpha =1)
8 plot(glm_lasso ,xvar="lambda",
col=colrs)
9 plot(glm_lasso ,col=colrs)
Before focusing on the lasso for binary variables, consider the standard ols-lasso...
But ﬁrst, to illustrate, let us get back on quantile ( 1) regression...

Optimization with linear problems (application to 1-regression)
Consider the case of the median. Consider a sample {y1, · · · , yn}.
To compute the median, solve min
µ
n
i=1
|yi − µ| which can be solved using linear
programming techniques (see e.g. simplex method)
More precisely, this problem is equivalent to min
µ,a,b
n
i=1
ai + bi with ai, bi ≥ 0 and
yi − µ = ai − bi, ∀i = 1, · · · , n.
Consider a sample obtained from a lognormal distribution
1 n = 101
2 set.seed (1)
3 y = rlnorm(n)
4 median(y)
5 [1] 1.077415

Here, one can use a standard optimization routine
1 md=Vectorize(function(m) sum(abs(y-m)))
2 optim(mean(y),md)
3 $par
4 [1] 1.077416
or a linear programing technique : use the matrix form, with 3n constraints, and 2n + 1
parameters,
1 library(lpSolve)
2 A1 = cbind(diag (2*n) ,0)
3 A2 = cbind(diag(n), -diag(n), 1)
4 r = lp("min", c(rep (1,2*n) ,0),
5 rbind(A1 , A2),c(rep(">=", 2*n), rep("=", n)), c(rep (0,2*n), y))
6 tail(r$solution ,1)
7 [1] 1.077415

More generally, consider here some quantile,
1 tau = .3
2 quantile(x,tau)
3 30%
4 0.6741586
The linear program is now min
µ,a,b
n
i=1
τai + (1 − τ)bi with ai, bi ≥ 0 and yi − µ = ai − bi,
∀i = 1, · · · , n.
1 A1 = cbind(diag (2*n) ,0)
2 A2 = cbind(diag(n), -diag(n), 1)
3 r = lp("min", c(rep(tau ,n),rep(1-tau ,n) ,0),
4 rbind(A1 , A2),c(rep(">=", 2*n), rep("=", n)), c(rep (0,2*n), y))
6 [1] 0.6741586

To illustrate numerical computations, use
1 base=read.table("http:// freakonometrics .free.fr/rent98_00. txt",header
=TRUE)
The linear program for the quantile regression is now
min
µ,a,b
n
i=1
τai + (1 − τ)bi
with ai, bi ≥ 0 and yi − [βτ
0 + βτ
1 xi] = ai − bi, ∀i = 1, · · · , n.
1 require(lpSolve)
2 tau = .3
3 n=nrow(base)
4 X = cbind( 1, base$area)
5 y = base$rent_euro
6 A1 = cbind(diag (2*n), 0,0)
7 A2 = cbind(diag(n), -diag(n), X)

1 r = lp("min",
2 c(rep(tau ,n), rep(1-tau ,n) ,0,0), rbind(A1 , A2),
3 c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y))
5 [1] 148.946864 3.289674
see
1 library(quantreg)
2 rq(rent_euro˜area , tau=tau , data=base)
3 Coefficients :
4 (Intercept) area
5 148.946864 3.289674

If we get back to the original Lasso approach, the goal was to solve
min
1
2n
n
i=1
[yi − (β0 + xi
T
β)]2
+ λ
j
|βj|
Since the intercept is not subject to the penalty. The first order condition is then
∂
∂β0
y − Xβ − β01 2
= (Xβ − y)T
1 + β0 1 2
= 0
i.e.
β0 =
1
n2
(Xβ − y)T
1
For other parameters, assuming tha KKT conditions are satisfied, we can check if 0 contains
the subdifferential at the minimum,
0 ∈ ∂
1
2
y − Xβ 2
+ λ β 1 =
1
2
y − Xβ 2
+ ∂(λ β 1 )

For the term on the left, we recognize
1
2
y − Xβ 2
= −XT
(y − Xβ) = −g
so that the previous equation can be writen
gk ∈ ∂(λ|βk|) =



{+λ} if βk > 0
{−λ} if βk < 0
(−λ, +λ) if βk = 0
i.e. if βk = 0, then gk = sign(βk) · λ.
We can split βj into a sum of its positive and negative parts by replacing βj with β+
j − β−
j
where β+
j , β−
j ≥ 0.

Then the Lasso problem becomes
− log L(β) + λ
j
(β+
j − β−
j )
with constraints β+
j − β−
j . Let α+
j , α−
j denote the Lagrange multipliers for β+
j , β−
j , respectively.
L(β) + λ
j
(β+
j − β−
j ) −
j
α+
j β+
j −
j
α−
j β−
j .
To satisfy the stationarity condition, we take the gradient of the Lagrangian with respect to β+
j
and set it to zero to obtain
L(β)j + λ − α+
j = 0
We do the same with respect to β−
j to obtain
− L(β)j + λ − α−
j = 0.

For conveniency, introduce the soft-thresholding function
S(z, γ) = sign(z) · (|z| − γ)+ =



z − γ if γ < |z| and z > 0
z + γ if γ < |z| and z < 0
0 if γ ≥ |z|
Noticing that the optimization problem
1
2
y − Xβ 2
2
+ λ β 1
can also be written
min
p
j=1
−βols
j · βj +
1
2
β2
j + λ|βj|
observe that βj,λ = S(βols
j , λ) which is a coordinate-wise update.

Now, if we consider a (slightly) more general problem, with weights in the ﬁrst part
min
1
2n
n
i=1
ωi[yi − (β0 + xi
T
β)]2
+ λ
j
|βj|
the coordinate-wise update becomes
βj,λ,ω = S(βω−ols
j , λ)
An alternative is to set
rj = y − β01 +
k=j
βkxk = y − y
(j)

so that the optimization problem can be written, equivalently
min
1
2n
p
j=1
[rj − βjxj]2
+ λ|βj|
hence
min
1
2n
p
j=1
β2
j xj − 2βjrj
T
xj + λ|βj|
and one gets
βj,λ =
1
xj
2
S(rj
T
xj, nλ)
or, if we develop
βj,λ =
1
i
x2
ij
S
i
xi,j[yi − y
(j)
i ], nλ

Again, if there are weights ω = (ωi), the coordinate-wise update becomes
βj,λ,ω =
1
i
ωix2
ij
S
i
ωixi,j[yi − y
(j)
i ], nλ
The code to compute this componentwise descent is
1 soft_ thresholding = function(x,a){
2 result = numeric(length(x))
3 result[which(x > a)] = a
4 result[which(x < -a)] <- x[which(x < -a)] + a
5 return(result)
and the code

1 lasso_coord_desc = function(X,y,beta ,lambda ,tol =1e-6, maxiter =1000){
2 beta = as.matrix(beta)
3 X = as.matrix(X)
4 omega = rep(1/length(y),length(y))
5 obj = numeric(length =( maxiter +1))
6 betalist = list(length(maxiter +1))
7 betalist [[1]] = beta
8 beta0list = numeric(length(maxiter +1))
9 beta0 = sum(y-X%*%beta)/(length(y))
10 beta0list [1] = beta0
11 for (j in 1: maxiter){
12 for (k in 1: length(beta)){
13 r = y - X[,-k]%*%beta[-k] - beta0*rep(1, length(y))
14 beta[k] = (1/sum(omega*X[,k]ˆ2))*soft_ thresholding (t(omega*r)
%*%X[,k],length(y)*lambda)
15 }

17 beta0list[j+1] = beta0
18 betalist [[j+1]] = beta
19 obj[j] = (1/2)*(1/length(y))*norm(omega*(y - X%*%beta -
20 beta0*rep(1, length(y))),’F’)ˆ2 + lambda*sum(abs(beta))
21 if (norm(rbind(beta0list[j],betalist [[j]]) - rbind(beta0 ,beta),
’F’) < tol) { break }
22 }
23 return(list(obj=obj [1:j],beta=beta ,intercept=beta0)) }
If we get back to our logistic problem, the log-likelihood is here
log L =
1
n
n
i=1
yi · (β0 + xi
T
β) − log[1 + exp(β0 + xi
T
β)]
which is a concave function of the parameters. Hence, one can use a quadratic approximation of
the log-likelihood - using Taylor expansion,
log L ≈ log L =
1
n
n
i=1
ωi · [zi − (β0 + xi
T
β)]2

where zi is the working response
zi = (β0 + xi
T
β) +
yi − pi
pi[1 − pi]
pi is the prediction
pi =
exp[β0 + xi
Tβ]
1 + exp[β0 + xi
Tβ]
and ωi are weights ωi = pi[1 − pi].
Thus, we obtain a penalized least-square problem.
1 lasso_coord_desc = function(X,y,beta ,lambda ,tol =1e-6, maxiter =1000){
2 beta = as.matrix(beta)
3 X = as.matrix(X)
4 obj = numeric(length =( maxiter +1))
5 betalist = list(length(maxiter +1))
6 betalist [[1]] = beta
8 p = exp(beta0*rep(1, length(y)) + X%*%beta)/(1+ exp(beta0*rep(1,
length(y)) + X%*%beta))
9 z = beta0*rep(1, length(y)) + X%*%beta + (y-p)/(p*(1-p))

10 omega = p*(1-p)/(sum((p*(1-p))))
11 beta0list = numeric(length(maxiter +1))
13 beta0list [1] = beta0
14 for (j in 1: maxiter){
15 for (k in 1: length(beta)){
16 r = z - X[,-k]%*%beta[-k] - beta0*rep(1, length(y))
17 beta[k] = (1/sum(omega*X[,k]ˆ2))*soft_ thresholding (t(omega*r)%
*%X[,k],length(y)*lambda)
18 }
20 beta0list[j+1] = beta0
21 betalist [[j+1]] = beta
22 obj[j] = (1/2)*(1/length(y))*norm(omega*(z - X%*%beta -
23 beta0*rep(1, length(y))),’F’)ˆ2 + lambda*sum(abs(beta))
24 p = exp(beta0*rep(1, length(y)) + X%*%beta)/(1+ exp(beta0*rep(1,
length(y)) + X%*%beta))
25 z = beta0*rep(1, length(y)) + X%*%beta + (y-p)/(p*(1-p))

26 omega = p*(1-p)/(sum((p*(1-p))))
27 if (norm(rbind(beta0list[j],betalist [[j]]) -
28 rbind(beta0 ,beta),’F’) < tol) { break }
29 }
30 return(list(obj=obj [1:j],beta=beta ,intercept=beta0)) }

4 df0 = df
5 df0$y = as.numeric(df$y) -1
6 plot_lambda = function(lambda){
7 m = apply(df0 ,2,mean)
8 s = apply(df0 ,2,sd)
9 for(j in 1:2) df0[,j] <- (df0[,j]-m[j])/s[j]
10 reg = glmnet(cbind(df0$x1 ,df0$x2), df0$y==1, alpha =1, lambda=lambda)
11 u = seq(0,1, length =101)
13 xt = (x-m[1])/s[1]
14 yt = (y-m[2])/s[2]
15 predict(reg ,newx=cbind(x1=xt ,x2=yt),type="response")}

4 reg = glmnet(cbind(df0$x1 ,df0$
x2), df0$y==1, alpha =1)
6 plot(reg ,xvar="lambda",col=c("
blue","red"),lwd =2)
7 abline(v=exp ( -2.8))
8 plot_lambda(exp ( -2.8))

4 reg = glmnet(cbind(df0$x1 ,df0$
x2), df0$y==1, alpha =1)
6 plot(reg ,xvar="lambda",col=c("
blue","red"),lwd =2)
7 abline(v=exp ( -2.1))
8 plot_lambda(exp ( -2.1))

Classiﬁcation : (Linear) Fisher’s Discrimination
Consider the follwing naive classiﬁcation rule
m (x) = argminy{P[Y = y|X = x]}
or
m (x) = argminy
P[X = x|Y = y]
P[X = x]
(where P[X = x] is the density in the continuous case).
In the case where y takes two values, that will be standard {0, 1} here, one can rewrite the later
as
m (x) =
1 if E(Y |X = x) >
1
2
0 otherwise

the set
DS = x, E(Y |X = x) =
1
2
is called the decision boundary.
Assume that X|Y = 0 ∼ N(µ0, Σ) and X|Y = 1 ∼ N(µ1, Σ), then explicit expressions can be
derived.
m (x) =
1 if r2
1 < r2
0 + 2log
P(Y = 1)
P(Y = 0)
+ log
|Σ0|
|Σ1|
0 otherwise
where r2
y is the Mahalanobis distance,
r2
y = [X − µy]T
Σ−1
y [X − µy]

Let δy be deﬁned as
δy(x) = −
1
2
log |Σy| −
1
2
[x − µy]T
Σ−1
y [x − µy] + log P(Y = y)
the decision boundary of this classiﬁer is
{x such that δ0(x) = δ1(x)}
which is quadratic in x. This is the quadratic discriminant analysis.
In Fisher (1936, The use of multile measurements in taxinomic problems , it was assumed that
Σ0 = Σ1
In that case, actually,
δy(x) = xT
Σ−1
µy −
1
2
µT
y Σ−1
µy + log P(Y = y)
and the decision frontier is now linear in x.

1 m0 = apply(myocarde[myocarde$PRONO =="0" ,1:7],2, mean)
2 m1 = apply(myocarde[myocarde$PRONO =="1" ,1:7],2, mean)
3 Sigma = var(myocarde [ ,1:7])
4 omega = solve(Sigma)%*%(m1 -m0)
5 omega
6 [,1]
7 FRCAR -0.012909708542
8 INCAR 1.088582058796
9 INSYS -0.019390084344
10 PRDIA -0.025817110020
11 PAPUL 0.020441287970
12 PVENT -0.038298291091
13 REPUL -0.001371677757
14 b = (t(m1)%*%solve(Sigma)%*%m1 -t(m0)%*%solve(Sigma)%*%m0)/2

1 x = c(.4 ,.55 ,.65 ,.9 ,.1 ,.35 ,.5 ,.15 ,.2 ,.85)
2 y = c(.85 ,.95 ,.8 ,.87 ,.5 ,.55 ,.5 ,.2 ,.1 ,.3)
3 z = c(1,1,1,1,1,0,0,1,0,0)
4 df = data.frame(x1=x,x2=y,y=as.factor(z))
5 m0 = apply(df[df$y=="0" ,1:2],2, mean)
6 m1 = apply(df[df$y=="1" ,1:2],2, mean)
7 Sigma = var(df [ ,1:2])
8 omega = solve(Sigma)%*%(m1 -m0)
9 omega
10 [,1]
11 x1 -2.640613174
12 x2 4.858705676

1 library(MASS)
2 fit_lda = lda(y ˜x1+x2 , data=df)
3 fit_lda
4
5 Coefficients of linear discriminants :
6 LD1
7 x1 -2.588389554
8 x2 4.762614663
and for the intercept
1 b = (t(m1)%*%solve(Sigma)%*%m1 -t(m0)%*%solve(
Sigma)%*%m0)/2

1 predlda = function(x,y) predict(fit_lda , data.
frame(x1=x,x2=y))$class ==1
2 vv=outer(vu ,vu ,predlda)
3 contour(vu ,vu ,vv ,add=TRUE ,lwd=2, levels = .5)
1 fit_qda = qda(y ˜x1+x2 , data=df)
2 plot(df$x1 ,df$x2 ,pch =19,
3 col=c("blue","red")[1+( df$y=="1")])
4 predqda=function(x,y) predict(fit_qda , data.
frame(x1=x,x2=y))$class ==1
5 vv=outer(vu ,vu ,predlda)
6 contour(vu ,vu ,vv ,add=TRUE ,lwd=2, levels = .5)

Note that there are strong links with the logistic regression.
Assume as previously that X|Y = 0 ∼ N(µ0, Σ) and X|Y = 1 ∼ N(µ1, Σ) then
log
P(Y = 1|X = x)
P(Y = 0|X = x)
is equal to
xT
Σ−1
[µy] −
1
2
[µ1 − µ0]T
Σ−1
[µ1 − µ0] + log
P(Y = 1)
P(Y = 0)
which is linear in x
log
P(Y = 1|X = x)
P(Y = 0|X = x)
= xT
β
Hence, when each groups have Gaussian distributions with identical variance matrix, then LDA
and the logistic regression lead to the same classiﬁcation rule.

Observe furthermore that the slope is proportional to Σ−1[µ1 − µ0] as stated in Fisher’s article.
But to obtain such a relationship, he observe that the ratio of between and within variances (in
the two groups) was
variance between
variance within
=
[ωµ1 − ωµ0]2
ωTΣ1ω + ωTΣ0ω
which is maximal when ω is proportional to Σ−1[µ1 − µ0].

Classiﬁcation : Neural Nets
The most interesting part is not ’neural’ but ’network’.

Networks have nodes, and edges (possibly connected) that connect nodes,
or maybe, to more speciﬁc some sort of ﬂow network,

In such a network, we usually have sources (here multiple) sources (here {s1, s2, s3}), on the
left, on a sink (here {t}), on the right. To continue with this metaphorical introduction,
information from the sources should reach the sink. An usually, sources are explanatory
variables, {x1, · · · , xp}, and the sink is our variable of interest y.
The output here is a binary variable y ∈ {0, 1}
Consider the sigmoid function f(x) =
1
1 + e−x
=
ex
ex + 1
which is actually the logistic function.
1 sigmoid = function(x) 1 / (1 + exp(-x))
This function f is called the activation function.

If y ∈ {−1, +1}, one can consider the hyperbolic tangent f(x) = tanh(x) =
(ex − e−x)
(ex + e−x)
or the
inverse tangent function f(x) = tan−1(x).
And as input for such function, we consider a weighted sum of incoming nodes.
So here
yi = f
p
j=1
ωjxj,i
We can also add a constant actually
yi = f ω0 +
p
j=1
ωjxj,i

1 weights_0 = lm(PRONO˜.,data=myocarde)$
coefficients
2 X = as.matrix(cbind (1, myocarde [ ,1:7]))
3 y_5_1 = sigmoid(X %*% weights_0)
4 library(ROCR)
5 pred = ROCR :: prediction(y_5_1,myocarde$PRONO)
6 perf = ROCR :: performance (pred ,"tpr", "fpr")
7 plot(perf ,col="blue",lwd =2)
8 reg = glm(PRONO˜.,data=myocarde ,family=
binomial(link = "logit"))
9 y_0 = predict(reg ,type="response")
10 pred0 = ROCR :: prediction(y_0,myocarde$PRONO)
11 perf0 = ROCR :: performance(pred0 ,"tpr", "fpr")
12 plot(perf0 ,add=TRUE ,col="red")

The error of that model can be computed used the quadratic loss function,
n
i=1
(yi − p(xi))2
or cross-entropy
n
i=1
(yi log p(xi) + [1 − yi] log[1 − p(xi)])
1 loss = function(weights){
2 mean( (myocarde$PRONO -sigmoid(X %*% weights))ˆ2) }
we want to solve
ω = argmin
1
n
n
i=1
yi, f(ω0 + xi
T
ω)

1 weights_1 = optim(weights_0,loss)$par
2 y_5_2 = sigmoid(X %*% weights_1)
6 plot(perf0 ,add=TRUE ,col="red")
Consider a single hidden layer, with 3 diﬀerent neurons, so
that
p(x) = f ω0 +
3
h=1
ωhfh ωh,0 +
p
j=1
ωh,jxj
or
p(x) = f ω0 +
3
h=1
ωhfh ωh,0 + xT
ωh

1 weights_1 = lm(PRONO˜1+ FRCAR+INCAR+INSYS+PAPUL
+PVENT ,data=myocarde)$ coefficients
2 X1 = as.matrix(cbind (1, myocarde[,c("FRCAR","
INCAR","INSYS","PAPUL","PVENT")]))
3 weights_2 = lm(PRONO˜1+ INSYS+PRDIA ,data=
myocarde)$ coefficients
4 X2 = as.matrix(cbind (1, myocarde[,c("INSYS","
PRDIA")]))
5 weights_3 = lm(PRONO˜1+ PAPUL+PVENT+REPUL ,data=
myocarde)$ coefficients
6 X3 = as.matrix(cbind (1, myocarde[,c("PAPUL","
PVENT","REPUL")]))
7 X = cbind(sigmoid(X1 %*% weights_1), sigmoid(
X2 %*% weights_2), sigmoid(X3 %*% weights_
3))
8 weights = c(1/3,1/3,1/3)
9 y_5_3 = sigmoid(X %*% weights)

If
p(x) = f ω0 +
3
h=1
ωhfh ωh,0 + xT
ωh
we want to solve, with back propagation,
ω = argmin
1
n
n
i=1
(yi, p(xi))
For practical issues, let us center all variables
1 center = function(z) (z-mean(z))/sd(z)
The loss function is here
1 loss = function(weights){
2 weights_1 = weights [0+(1:7)]
5 weights_ = weights [21+1:4]
6 X1=X2=X3=as.matrix(myocarde [ ,1:7])
7 Z1 = center(X1 %*% weights_1)

10 X = cbind (1, sigmoid(Z1), sigmoid(Z2), sigmoid(Z3))
11 mean( (myocarde$PRONO -sigmoid(X %*% weights_))ˆ2)}

1 pca = princomp(myocarde [ ,1:7])
2 W = get_pca_var(pca)$contrib
3 weights_0 = c(W[,1],W[,2],W[,3],c(-1,rep (1 ,3)/
3))
4 weights_opt = optim(weights_0,loss)$par
5 weights_1 = weights_opt [0+(1:7)]
8 weights_ = weights_opt [21+1:4]
9 X1=X2=X3=as.matrix(myocarde [ ,1:7])
13 X = cbind (1, sigmoid(Z1), sigmoid(Z2), sigmoid(
Z3))
14 y_5_4 = sigmoid(X %*% weights_)
18 plot(perf ,add=TRUE ,col="red")

1 library(nnet)
2 myocarde_minmax = myocarde
3 minmax = function(z) (z-min(z))/(max(z)-min(z)
)
4 for(j in 1:7) myocarde_minmax[,j] = minmax(
myocarde_minmax[,j])
5 model_nnet = nnet(PRONO˜.,data=myocarde_minmax
,size =3)
6 summary(model_nnet)
7 library(devtools)
8 source_url(’https://gist. githubusercontent .com
/fawda123/7471137/raw/466
c1474d0a505ff044412703516c34f1a4684a5 /nnet
_plot_update.r’)
9 plot.nnet(model_nnet)
10 library(neuralnet)
11 model_nnet = neuralnet(formula(glm(PRONO˜.,
data=myocarde_minmax)),
12 myocarde_minmax ,hidden =3, act.fct = sigmoid)
13 plot(model_nnet)

Of course, one can consider a multiple layer network,
p(x) = f ω0 +
3
h=1
ωhfh ωh,0 + zh
T
ωh
where
zh = f ωh,0 +
kh
j=1
ωh,jfh,j ωh,j,0 + xT
ωh,j
1 library(neuralnet)
2 model_nnet = neuralnet(formula(glm(PRONO˜.,
data=myocarde_minmax)),
3 myocarde_minmax ,hidden =3, act.fct = sigmoid)
4 plot(model_nnet)

SVM : Support Vector Machine
Linearly Separable sample [econometric notations]
Data (y1, x1), · · · , (yn, xn) - with y ∈ {0, 1} - are
linearly separable if there are (β0, β) such that
- yi = 1 if β0 + x β > 0
- yi = 0 if β0 + x β < 0
Linearly Separable sample [ML notations]
Data (y1, x1), · · · , (yn, xn) - with y ∈ {−1, +1} -
are linearly separable if there are (b, ω) such that
- yi = +1 if b + x, ω > 0
- yi = −1 if b + x, ω < 0
or equivalently yi · (b + x, ω ) > 0, ∀i.

yi · (b + x, ω ) = 0 is an hyperplane (in Rp)
orthogonal with ω
Use m(x) = 1b+ x,ω ≥0 − 1b+ x,ω <0 as classiﬁer
Problem : equation (i.e. (b, ω)) is not unique !
Canonical form : min
i=1,··· ,n
|b + xi, ω | = 1
Problem : solution here is not unique !
Idea : use the widest (safety) margin γ
Vapnik & Lerner (1963, Pattern recognition using generalized
portrait method) or Cover (1965, Geometrical and statistical
properties of systems of linear inequalities with applications in
pattern recognition)

Consider two points, ω−1 and ω+1
γ =
1
2
ω, ω+1 − ω−1
ω
It is minimal when
b + xi, ω−1 = −1 and
b + xi, ω+1 = +1, and therefore
γ =
1
ω
Optimization problem becomes
min
(b ω)
1
2
ω 2
2
s.t. yi · (b + x, ω ) > 0, ∀i.
convex optimization problem with linear constraints

Consider the following problem : min
u∈Rp
h(u) s.t. gi(u) ≥ 0 ∀i = 1, · · · , n
where h is quadratic and gi’s are linear.
Lagrangian is L : Rp × Rn → R deﬁned as L(u, α) = h(u) −
n
i=1
αigi(u)
where α are dual variables, and the dual function is
Λ : α → L(uα, α) = min{L(u, α)} where uα = argmin L(u, α)
One can solve the dual problem, max{Λ(α)} s.t. α ≥ 0. Solution is u = uα .
Si gi(uα ) > 0, then necessarily αi = 0 (see Karush-Kuhn-Tucker (KKT) condition,
αi · gi(uα ) = 0)

Here, L(b, ω, α) =
1
2
ω 2
−
n
i=1
αi · yi · (b + x, ω ) − 1
From the ﬁrst order conditions,
∂L(b, ω, α)
∂ω
= ω −
n
i=1
αi · yixi = 0, i.e.ωËĘ =
n
i=1
αi yixi
∂L(b, ω, α)
∂b
= −
n
i=1
αi · yi = 0, i.e.
n
i=1
αi · yi
and
Λ(α) =
n
i=1
αi −
1
2
n
i,j=1
αiαjyiyj αi, αj

min
α
1
2
α Qα − 1 α s.t.
αi ≥ 0, ∀i
y 1 = 0
where Q = [Qi,j] and Qi,j = yiyj xi, xj , and then
ω =
n
i=1
αi yixi and b = −
1
2
min
i:yi=+1
{ xi, ω } + min
i:yi=−1
{ xi, ω }
Points xi such that αi > 0 are called support
yi · (b + xi, ω ) = 1
Use m (x) = 1b + x,ω ≥0 − 1b + x,ω <0 as classiﬁer
Observe that γ =
n
i=1
α 2
i
−1/2

Optimization problem was
min
(b,ω)
1
2
ω 2
2
s.t. yi · (b + x, ω ) > 0, ∀i,
which became
min
(b,ω)
1
2
ω 2
2
+
n
i=1
αi · 1 − yi · (b + x, ω ) ,
or
min
(b,ω)
1
2
ω 2
2
+ penalty ,

Consider here the more general case
where the space is not linearly separable
( ω, xi + b)yi ≥ 1
becomes
( ω, xi + b)yi ≥ 1−ξi
for some slack variables ξi’s.
and penalize large slack variables ξi, i.e. solve (for some coste C)
min
ω,b
1
2
ω β + C
n
i=1
ξi
subject to ∀i, ξi ≥ 0 and (xi ω + b)yi ≥ 1 − ξi.
This is the soft-margin extension, see e1071::svm() or kernlab::ksvm()

The dual optimization problem is now
min
α
1
2
α Qα − 1 α s.t.
0 ≤ αi≤ C, ∀i
y 1 = 0
where Q = [Qi,j] and Qi,j = yiyj xi, xj , and then
ω =
n
i=1
αi yixi and b = −
1
2
min
i:yi=+1
{ xi, ω } + min
i:yi=−1
{ xi, ω }
Note further that the (primal) optimization problem can be written
min
(b,ω)
1
2
ω 2
2
+
n
i=1
1 − yi · (b + x, ω )
+
,
where (1 − z)+ is a convex upper bound for empirical error 1z≤0

One can also consider the kernel trick : xi xj is replace by ϕ(xi) ϕ(xj) for some mapping ϕ,
K(xi, xj) = ϕ(xi) ϕ(xj)
For instance
K(a, b) = (a b)3
= ϕ(a) ϕ(b)
where
ϕ(a1, a2) = (a3
1 ,
√
3a2
1a2 ,
√
3a1a2
2 , a3
2)
Consider polynomial kernels
K(a, b) = (1 + a b)p
or a Gaussian kernel
K(a, b) = exp(−(a − b) (a − b))
Solve
max
αi≥0
n
i=1
αi −
1
2
n
i,j=1
yiyjαiαjK(xi, xj)

Linear kernel

Polynomial kernel (degree 2)

Radial kernel

SVM : Support Vector Machine - Radial Kernel, impact of the cost C

SVM : Support Vector Machine - Radial Kernel, tuning parameter γ

The radial kernel is formed by taking an inﬁnite sum over polynomial kernels...
K(x, y) = exp − γ x − y 2
= ψ(x), ψ(y)
where ψ is some Rn → R∞ function, since
K(x, y) = exp − γ x − y 2
= exp(−γ x 2
− γ y 2
)
=constant
· exp 2γ x, y
i.e.
K(x, y) ∝ exp 2γ x, y =
∞
k=0
2γ
x, y k
k!
=
∞
k=0
2γKk(x, y)
where Kk is the polynomial kernel of degree k.
If K = K1 + K2 with ψj : Rn → Rdj then ψ : Rn → Rd with d ∼ d1 + d2

A kernel is a measure of similarity between vectors.
The smaller the value of γ the narrower the vectors should be to have a small measure
Is there a probabilistic interpretation ?
Platt (2000, Probabilities for SVM) suggested to use a logistic function over the SVM scores,
p(x) =
exp[b + x, ω ]
1 + exp[b + x, ω ]

Let C denote the cost of mis-classiﬁcation. The optimization problem becomes
min
w∈Rd,b∈R,ξ∈Rn
1
2
ω 2
+ C
n
i=1
ξi
subject to yi · (ωTxi + b) ≥ 1 − ξi, with ξi ≥ 0, ∀i = 1, · · · , n
1 n = length(myocarde[,"PRONO"])
2 myocarde0 = myocarde
3 myocarde0$PRONO = myocarde$PRONO*2-1
4 C = .5
5 f = function(param){
6 w = param [1:7]
7 b = param [8]
8 xi = param [8+1: nrow(myocarde)]
9 .5*sum(wˆ2) + C*sum(xi)}
10 grad_f = function(param){
11 w = param [1:7]
12 b = param [8]

13 xi = param [8+1: nrow(myocarde)]
14 c(2*w,0,rep(C,length(xi)))}
(linear) constraints are written as Uθ − c ≥ 0
1 U = rbind(cbind(myocarde0[,"PRONO"]*as.matrix(myocarde [ ,1:7]) ,diag(n)
,myocarde0[,"PRONO"]),
2 cbind(matrix (0,n,7),diag(n,n),matrix (0,n,1)))
3 C = c(rep(1,n),rep(0,n))
then use
1 constrOptim (theta=p_init , f, grad_f, ui = U,ci = C)
One can recognize in the separable case, but also in the non-separable case, a classic quadratic
program
min
z∈Rd
1
2
zT
Dz − dz
subject to Az ≥ b.
1 library(quadprog)
2 eps = 5e-4

3 y = myocarde[,"PRONO&quot ;]*2-1
5 n = length(y)
6 D = diag(n+7+1)
7 diag(D)[8+0:n] = 0
8 d = matrix(c(rep (0 ,7) ,0,rep(C,n)), nrow=n+7+1)
9 A = Ui
10 b = Ci

Classiﬁcation : Support Vector Machine (SVM) - Primal Problem
1 sol = solve.QP(D+eps*diag(n+7+1) , d, t(A), b, meq =1)
2 qpsol = sol$solution
3 (omega = qpsol [1:7])
4 [1] -0.106642005446 -0.002026198103 -0.022513312261 -0.018958578746
-0.023105767847 -0.018958578746 -1.080638988521
5 (b = qpsol[n+7+1])
6 [1] 997.6289927
Given an observation x, the prediction is y = sign[ωTx + b]
1 y_pred = 2*((as.matrix(myocarde0 [ ,1:7])%*%omega+b) >0) -1

Classiﬁcation : Support Vector Machine (SVM) - Dual Problem
The Lagrangian of the separable problem could be written introducing Lagrange multipliers
α ∈ Rn, α ≥ 0 as
L(ω, b, α) =
1
2
ω 2
−
n
i=1
αi yi(ωT
xi + b) − 1
Somehow, αi represents the inﬂuence of the observation (yi, xi). Consider the Dual Problem,
with G = [Gij] and Gij = yiyjxj
Txi.
min
α∈Rn
1
2
αT
Gα − 1T
α
subject to yTα = 0.
The Lagrangian of the non-separable problem could be written introducing Lagrange
multipliers α, β ∈ Rn, α, β ≥ 0,

and deﬁne the Lagrangian L(ω, b, ξ, α, β) as
1
2
ω 2
+ C
n
i=1
ξi −
n
i=1
αi yi(ωT
xi + b) − 1 + ξi −
n
i=1
βiξi
Somehow, αi represents the inﬂuence of the observation (yi, xi). The Dual Problem become
with G = [Gij] and Gij = yiyjxj
Txi.
min
α∈Rn
1
2
αT
Gα − 1T
α
subject to yTα = 0, α ≥ 0 and α ≤ C. As previously, one can also use quadratic programming

1 library(quadprog)
2 eps = 5e-4
3 y = myocarde[,"PRONO"]*2-1
5 n = length(y)
6 Q = sapply (1:n, function(i) y[i]*t(X)[,i])
7 D = t(Q)%*%Q
8 d = matrix (1, nrow=n)
9 A = rbind(y,diag(n),-diag(n))
10 C = .5
11 b = c(0,rep(0,n),rep(-C,n))
12 sol = solve.QP(D+eps*diag(n), d, t(A), b, meq=1, factorized=FALSE)
13 qpsol = sol$solution

The two problems are directed in the sense that for all xωTx + b =
n
i=1
αiyi(xTxi) + b
To recover the solution of the primal problem, ω =
n
i=1
αiyixi thus
1 omega = apply(qpsol*y*X,2,sum)
2 omega
3 1 FRCAR INCAR INSYS
4 0.000000000000000243907 0.0550138658 -0.092016323 0.3609571899422
5 PRDIA PAPUL PVENT REPUL
6 -0.109401796528869235669 -0.0485213403 -0.066005864 0.0010093656567
More generally
1 svm.fit = function(X, y, C=NULL) {
2 n.samples = nrow(X)
3 n.features = ncol(X)
4 K = matrix(rep(0, n.samples*n.samples), nrow=n.samples)
5 for (i in 1:n.samples){
6 for (j in 1:n.samples){
7 K[i,j] = X[i,] %*% X[j,] }}

8 Dmat = outer(y,y) * K
9 Dmat = as.matrix(nearPD(Dmat)$mat)
10 dvec = rep(1, n.samples)
11 Amat = rbind(y, diag(n.samples), -1*diag(n.samples))
12 bvec = c(0, rep(0, n.samples), rep(-C, n.samples))
13 res = solve.QP(Dmat ,dvec ,t(Amat),bvec=bvec , meq =1)
14 a = res$solution
15 bomega = apply(a*y*X,2,sum)
16 return(bomega)
17 }

To relate this duality optimization problem to OLS, recall that y = xTω + ε, so that y = xTω,
where ω = [XTX]−1XTy.
But one can also write y = xTω =
n
i=1
αi · xTxi where α = X[XTX]−1ω, or conversely,
ω = XTα.

Classiﬁcation : Support Vector Machine (SVM) - Kernel Approach
A positive kernel on X is a function K : X × X → R symmetric, and such that for any n,
∀α1, · · · , αn and ∀x1, · · · , xn,
n
i=1
n
j=1
αiαjk(xi, xj) ≥ 0
For example, the linear kernel is k(xi, xj) = xi
Txj. ThatâĂŹs what weâĂŹve been using here,
so far. One can also deﬁne the product kernel k(xi, xj) = κ(xi) · κ(xj) where κ is some
function X → R.

Classiﬁcation : Support Vector Machine (SVM) - Kernel Approach
Finally, the Gaussian kernel is k(xi, xj) = exp[− xi − xj
2].
Since it is a function of xi − xj , it is also called a radial kernel.
1 linear.kernel = function(x1 , x2) {
2 return (x1%*%x2)
3 }
4 svm.fit = function(X, y, FUN=linear.kernel , C=NULL) {
5 n.samples = nrow(X)

6 n.features = ncol(X)
7 K = matrix(rep(0, n.samples*n.samples), nrow=n.samples)
8 for (i in 1:n.samples){
9 for (j in 1:n.samples){
10 K[i,j] = FUN(X[i,], X[j,])
11 }
12 }
13 Dmat = outer(y,y) * K
14 Dmat = as.matrix(nearPD(Dmat)$mat)
15 dvec = rep(1, n.samples)
16 Amat = rbind(y, diag(n.samples), -1*diag(n.samples))
17 bvec = c(0, rep(0, n.samples), rep(-C, n.samples))
18 res = solve.QP(Dmat ,dvec ,t(Amat),bvec=bvec , meq =1)
19 a = res$solution
20 bomega = apply(a*y*X,2,sum)
21 return(bomega)
22 }

1 SVM2 = ksvm(y ˜ x1 + x2 , data = df0 , C=2,
kernel = "vanilladot " , prob.model=TRUE ,
type="C-svc")
2 pred_SVM2 = function(x,y){
3 return(predict(SVM2 ,newdata=data.frame(x1=x,
x2=y), type=" probabilities ")[ ,2])}
4 plot(df$x1 ,df$x2 ,pch=c(1 ,19) [1+( df$y=="1")],
5 cex =1.5 , xlab="",
6 ylab="",xlim=c(0 ,1),ylim=c(0 ,1))
7 vu = seq (-.1,1.1, length =251)
8 vv = outer(vu ,vu ,function(x,y) pred_SVM2(x,y
))
9 contour(vu ,vu ,vv ,add=TRUE ,lwd=2, levels = .5,
col="red")

1 SVM3 = ksvm(y ˜ x1 + x2 , data = df0 , C=2,
kernel = "rbfdot" , prob.model=TRUE ,
type="C-svc")
2 pred_SVM3 = function(x,y){
3 return(predict(SVM3 ,newdata=data.frame(x1=x,
x2=y), type=" probabilities ")[ ,2])}
4 plot(df$x1 ,df$x2 ,pch=c(1 ,19) [1+( df$y=="1")],
5 cex =1.5 , xlab="",
6 ylab="",xlim=c(0 ,1),ylim=c(0 ,1))
7 vu = seq (-.1,1.1, length =251)
8 vv = outer(vu ,vu ,function(x,y) pred_SVM2(x,y
))
9 contour(vu ,vu ,vv ,add=TRUE ,lwd=2, levels = .5,
col="red")

Classification : Classification Trees
1 library(rpart)
2 cart = rpart(PRONO˜.,data=
myocarde)
3 library(rpart.plot)
4 prp(cart ,type=2, extra =1)
A (binary) split is based on one specific variable âĂŞ say xj âĂŞ and a cutoff, say s. Then,
there are two options:
• either xi,j ≤ s, then observation i goes on the left, in IL
• or xi,j > s, then observation i goes on the right, in IR
Thus, I = IL ∪ IR.

Gini for node I is deﬁned as
G(I) = −
y∈{0,1}
py(1 − py)
where py is the proportion of individuals in the leaf of type y,
G(I) = −
y∈{0,1}
ny,I
nI
1 −
ny,I
nI
1 gini = function(y,classe){
2 T. = table(y,classe)
3 nx = apply(T,2,sum)
4 n. = sum(T)
5 pxy = T/matrix(rep(nx ,each =2) ,nrow =2)
6 omega = matrix(rep(nx ,each =2) ,nrow =2)/n
7 g. = -sum(omega*pxy*(1-pxy))
8 return(g)}

1 -2*mean(myocarde$PRONO)*(1-mean(myocarde$PRONO))
2 [1] -0.4832375
3 gini(y=myocarde$PRONO ,classe=myocarde$PRONO <Inf)
4 [1] -0.4832375
5 gini(y=myocarde$PRONO ,classe=myocarde [ ,1] <=100)
6 [1] -0.4640415

if we split, deﬁne index
G(IL, IR) = −
x∈{L,R}
nx
nIx
nI
y∈{0,1}
ny,Ix
nIx
1 −
ny,Ix
nIx
the entropic measure is
E(I) = −
y∈{0,1}
ny,I
nI
log
ny,I
nI
1 entropy = function(y,classe){
2 T = table(y,classe)
3 nx = apply(T,2,sum)
4 pxy = T/matrix(rep(nx ,each =2) ,nrow =2)
5 omega = matrix(rep(nx ,each =2) ,nrow =2)/sum(T)
6 g = sum(omega*pxy*log(pxy))
7 return(g)}

1 mat_gini = mat_v=matrix(NA ,7 ,101)
2 for(v in 1:7){
3 variable=myocarde[,v]
4 v_seuil=seq(quantile(myocarde[,v],
5 6/length(myocarde[,v])),
6 quantile(myocarde[,v],1-6/length(
7 myocarde[,v])),length =101)
8 mat_v[v,]=v_seuil
9 for(i in 1:101){
10 CLASSE=variable <=v_seuil[i]
11 mat_gini[v,i]=
12 gini(y=myocarde$PRONO ,classe=CLASSE)}}
13 -(gini(y=myocarde$PRONO ,classe =( myocarde
[ ,3] <19))-
14 gini(y=myocarde$PRONO ,classe =( myocarde [,3]<
Inf)))/
15 gini(y=myocarde$PRONO ,classe =( myocarde [,3]<
Inf))
16 [1] 0.5862131

1 idx = which(myocarde$INSYS <19)
2 mat_gini = mat_v = matrix(NA ,7 ,101)
3 for(v in 1:7){
4 variable = myocarde[idx ,v]
5 v_seuil = seq(quantile(myocarde[idx ,v],
6 7/length(myocarde[idx ,v])),
7 quantile(myocarde[idx ,v],1-7/length(
8 myocarde[idx ,v])), length =101)
9 mat_v[v,] = v_seuil
10 for(i in 1:101){
11 CLASSE = variable <=v_seuil[i]
12 mat_gini[v,i]=
13 gini(y=myocarde$PRONO[idx],classe=
CLASSE)}}
15 for(v in 2:7){
16 plot(mat_v[v,],mat_gini[v ,])
17 }

1 idx = which(myocarde$INSYS >=19)
2 mat_gini = mat_v = matrix(NA ,7 ,101)
3 for(v in 1:7){
4 variable=myocarde[idx ,v]
5 v_seuil=seq(quantile(myocarde[idx ,v],
6 6/length(myocarde[idx ,v])),
7 quantile(myocarde[idx ,v],1-6/length(
8 myocarde[idx ,v])), length =101)
9 mat_v[v,]=v_seuil
10 for(i in 1:101){
11 CLASSE=variable <=v_seuil[i]
12 mat_gini[v,i]=
13 gini(y=myocarde$PRONO[idx],
14 classe=CLASSE)}}
16 for(v in 2:7){
17 plot(mat_v[v,],mat_gini[v ,])
18 }

Classiﬁcation : Variable Importance, on Trees
1 cart = rpart(PRONO˜., myocarde)
2 (split = summary(cart)$splits)
3 count ncat improve index adj
4 INSYS 71 -1 0.58621312 18.850 0.0000000
5 REPUL 71 1 0.55440034 1094.500 0.0000000
6 INCAR 71 -1 0.54257020 1.690 0.0000000
7 PRDIA 71 1 0.27284114 17.000 0.0000000
8 PAPUL 71 1 0.20466714 23.250 0.0000000
1 REPUL 27 1 0.18181818 1585.000 0.0000000
2 PVENT 27 -1 0.10803571 14.500 0.0000000
3 PRDIA 27 1 0.10803571 18.500 0.0000000
4 PAPUL 27 1 0.10803571 22.500 0.0000000
5 INCAR 27 1 0.04705882 1.195 0.0000000
1 cart$variable.importance
2 INSYS REPUL INCAR PAPUL PRDIA FRCAR PVENT
3 10.364984 10.05108 8.21212 3.24415 2.82761 1.86230 0.33737

Classiﬁcation : Bagging Trees
Random forests are obtained using
1 L_tree = list ()
2 for(s in 1:1000){
3 idx = sample (1:n, size=n, replace=TRUE)
4 L_tree [[s]] = rpart(as.factor(PRONO)˜.,
myocarde[idx ,])
5 }
Actually, random forest are slightly diﬀerent : at each node,
variables are randomly selected.

Boosting & Adaboost
Classiﬁcation problem, yi ∈ {•, •}, consider a model at stage
k − 1,
if mk−1(xi) = yi, increase the weight given to observation i
Boosting : weak learner
A weak model is a model slightly better than a
pure random one (head/tails)
q
q
q
q
q
q
q
q
q
q
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
0.5
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
0.0 0.2 0.4 0.6 0.8 1.0
0.00.20.40.60.81.0
0.5q
q
q
q
q
q
q
q
q

Boosting & Adaboost
Adaboost Algorithm
1. Set weights ωi = 1/n, i = 1, · · · , n
2 . For k = 1, · · ·
(i) ﬁt model on (yi, xi) with weights ωi, get hk(x)
(ii) compute the error rate εk =
n
i=1
˜ωi1yi=hk(xi)
(iii) compute αk = log
1 − εk
εk
(iv) reevaluate the weights ωi = ωi · eαk1yi=hk(xi)
3. The ﬁnal model is hκ(x) =
κ
k≥1
αkhk(x)
The error rate should not be too small (εk ≤ 50%) to insure αk > 0

Boosting & Adaboost
The general problem in machine learning is to ﬁnd m (·) = argmin
m∈M
E (Y, g(X)
Use loss (y, g(x)) = 1y=g(x.
Empirical version is mn(·) = argmin
m∈M
1
n
n
i=1
(yi, g(xi) = argmin
m∈M
1
n
n
i=1
1yi=g(xi)
Complicated problem : use a convex version of the loss function
(y, g(x) = exp[−y · g(x)]
From Hastie et al. (2009), with the adaboost algorithm,
hκ(·) = hκ−1(·) + ακhκ(x) = hκ−1(·) + 2β H (·)
where (β , H (·)) = argmin
(β,H)∈(R,M)
n
i=1
exp − yi · (hκ−1(xi) + βH(xi)

Boosting & Adaboost
From Freund & Shapire (1999), empirical error of hκ(·) satisﬁes
1
n
n
i=1
1yi=hκ(xi) ≤ exp −2
κ
k=1
(εk − 0.5)2
(when weak learners are better than random classiﬁcation, empirical error tends to 0 -
exponentially fast)
−3 −2 −1 0 1 2 3
0.00.20.40.60.81.0
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
qq
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
qq
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
qq
q
q
q
q
q
q
q
q
q
−3 −2 −1 0 1 2 3
0.00.20.40.60.81.0
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
qq
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
q
qq
q
q
q
q
0 2000 4000 6000 8000
0.00.10.20.30.40.5
Number of iterations
Error
training sample
validation sample

Gradient Boosting
Newton-Raphson to minimize a strictly convex function g : R → R
At minimum, g (x ) = 0, so consider ﬁrst order approximation g (x + h) ≈ g (x) + h · g (x)
Consider sequence xk = xk−1 − αg (xk−1) where α = [g (xk−1)]−1
One can consider a functional version of that technique, ∀i = 1, · · · , n,
gk(xi) = gk−1(xi) − α
∂ (yi, g(xi))
∂g(xi) g(xi)=gk−1(xi)
This provides a sequence of function gk at points xi.
To get values at any point x use regression i’s on xi’s, εi = −
∂ (yi, g))
∂g g=gk−1(xi)
If α = 1 and (y, g) = exp[−yg], we have (almost) adaboost

Gradient Boosting
Gradient Boosting Algorithm
1. Start with a constant model, h0(x) = argmin
c∈R
1
n
n
i=1
(yi, c)
and a regularization parameter α ∈ (0, 1)
2 . For k = 1, · · ·
(i) compute εi = −
∂ (yi, g))
∂g g=gk−1(xi)
(ii) ﬁt the (weak) model on sample (εi, xi) and let Hk denote
the mode
(iii) update the model hk(·) = hk−1(·) + αHk(·)
3. The ﬁnal model is hκ(x)
The choice of α is (somehow) not important : use α ∼ 10%

Gradient Boosting
The logiboost model is obtained when y ∈ {0, 1} and loss function is
(y, m) = log[1 + exp(−2(2y − 1)m)]
Boosting (learning from the mistakes)
Sequential Learning
mk(·) = mk−1(·) + α · argmin
h∈H



n
i=1
yi − mk−1(xi)
εi
, h(xi)



Hence, learning is sequential, as opposed to bagging...

1 n_iter = 100
2 y = (myocarde[,"PRONO"]==1)*2-1
3 x = myocarde [ ,1:7]
4 error = rep(0,n_iter)
5 f = rep(0, length(y))
6 w = rep(1, length(y)) #
7 alpha = 1
8 library(rpart)
9 for(i in 1:n_iter){
10 w = exp(-alpha*y*f) *w
11 w = w/sum(w)
12 rfit = rpart(y˜., x, w, method="class")
13 g = -1 + 2*(predict(rfit ,x)[ ,2] >.5)
14 e = sum(w*(y*g<0))
15 alpha = .5*log ( (1-e) / e )
16 alpha = 0.1*alpha
17 f = f + alpha*g
18 error[i] = mean (1*f*y <0)
19 }
20 plot(seq(1,n_iter),error)

Classiﬁcation : Boosting (adabost)
It is necessary to use a training / validation dataset to avoid overﬁt...
1 set.seed (123)
2 id_train = sample (1: nrow(myocarde), size =45, replace=FALSE)
3 train_myocarde = myocarde[id_train ,]
4 test_myocarde = myocarde[-id_train ,]

1 for(i in 1:n_iter){
2 w_train = w_train*exp(-alpha*y_train*f_
train)
3 w_train = w_train/sum(w_train)
4 rfit = rpart(y_train˜., x_train , w_train ,
method="class")
5 g_train = -1 + 2*(predict(rfit ,x_train)
[ ,2] >.5)
6 g_test = -1 + 2*(predict(rfit ,x_test)
[ ,2] >.5)
7 e_train = sum(w_train*(y_train*g_train <0))
8 alpha = .5*log ( (1-e_train) / e_train )
9 alpha = 0.1*alpha
10 f_train = f_train + alpha*g_train
11 f_test = f_test + alpha*g_test
12 train_error[i] = mean (1*f_train*y_train <0)
13 test_error[i] = mean (1*f_test*y_test <0)}

Classiﬁcation : Boosting (adabost)
1 library(gbm)
2 gbmWithCrossValidation = gbm(PRONO ˜ .,
distribution = "bernoulli",
3 data = myocarde ,n.trees = 2000 , shrinkage =
.01,cv.folds = 5,n.cores = 1)
4 bestTreeForPrediction = gbm.perf(
gbmWithCrossValidation )

Gradient Boosting
We have used (y, g(x) = exp[−y · m(x)] instead of 1y=m(x.
Misclassification error is (upper) bounded by the exponential loss
1
n
n
i=1
1yi·m(xi
≤
1
n
n
i=1
exp[−yi · m(xi]
Here m(x) is a linear combination of weak classifier, m(x) =
κ
j=1
αjhj(x).
Let M = [Mi,j] where Mi,j = yi · hj(xi) ∈ {−1, +1}, i.e. Mi,j = 1 whenever (weak) classifier j
correctly classifies individual i.
yi · m(xi) =
κ
j=1
αjyihj(xi) = Mα
i
thus, R(α) =
1
n
n
i=1
exp[−yi · m(xi)] =
1
n
n
i=1
exp − (Mα)i

Gradient Boosting
One can use coordinate descent, in direction j in which the directional derivative is the steepest,
j ∈ argmin −
∂R(α + aej)
∂a a=0
where the objective can be written
−
∂
∂a
1
n
n
i=1
exp − (Mα)i − a(Mej)i
a=0
=
1
n
n
i=1
Mij exp − (Mα)i
Then
j ∈ argmin (d M)j where di =
exp[−(Mα)i]
i
exp[−(Mα)i]
Then do a line-search to see how far we should me in that direction. The derivative is null if
−
∂R(α + aej)
∂a
= 0 i.e. a =
1
2
log
d+
=
1
2
log
1 − d−
d−
where
d− =
i:Mi,j =−1
di and d+ =
i:Mi,j =+1
di

Gradient Boosting
The algorithm is here
Coordinate Descent Algorithm
1. di = 1/n for i = 1, · · · , n and α = 0
2 . For k = 1, · · ·
(i) ﬁnd optimal direction j ∈ argmin (d M)j
(ii) compute − =
i:Mi,j =−1
di and ak =
1
2
log
1 − d−
d−
(iii) set α = α + akej and di = exp[−(Mα)i]
i
exp[−(Mα)i]

Gradient Boosting
very close to Adaboost : αj is the sum of ak where direction j was considered,
αj =
κ
k=1
ak1j (k)=j
Thus
m (x) =
κ
k=1
αjhj(x) =
κ
k=1
akhj (k)(x)
With Adaboost, we go in the same direction, with the same intensity : Adaboost is equivalent
to minimizing the exponential loss by coordinate descent.
Thus, we seek m (·) = argmin E(Y,X)∼F
exp (−Y · m(X))
which is minimized at m (x) =
1
2
log
P[Y = +1|X = x]
P[Y = −1|X = x]
(very close to the logistic regression)

Hands-On Algorithms for Predictive Modeling

Recommended

More Related Content

What's hot (20)

Similar to Hands-On Algorithms for Predictive Modeling (20)

More from Arthur Charpentier (12)

Recently uploaded (20)

Hands-On Algorithms for Predictive Modeling