François presented our joint paper “A Sequentially Fair Mechanism for Multiple Sensitive Attributes“, written with Philipp Ratz, at the L^P seminar(ISFA-CNAM-HEC Lausanne) today
Thanks François ! (more to come soon)
François presented our joint paper “A Sequentially Fair Mechanism for Multiple Sensitive Attributes“, written with Philipp Ratz, at the L^P seminar(ISFA-CNAM-HEC Lausanne) today
Thanks François ! (more to come soon)
[an updated version is now online here]
After my series of post on classification algorithms, it’s time to get back to R codes, this time for quantile regression. Yes, I still want to get a better understanding of optimization routines, in R. Before looking at the quantile regression, let us compute the median, or the quantile, from a sample.
Consider a sample \{y_1,\cdots,y_n\}. To compute the median, solve\min_\mu \left\lbrace\sum_{i=1}^n|y_i-\mu|\right\rbracewhich can be solved using linear programming techniques. More precisely, this problem is equivalent to\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^na_i+b_i\right\rbracewith a_i,b_i\geq 0 and y_i-\mu=a_i-b_i, \forall i=1,\cdots,n.
To illustrate, consider a sample from a lognormal distribution,
n = 101 set.seed(1) y = rlnorm(n) median(y) [1] 1.077415 |
For the optimization problem, use the matrix form, with 3n constraints, and 2n+1 parameters,
library(lpSolve) A1 = cbind(diag(2*n),0) A2 = cbind(diag(n), -diag(n), 1) r = lp("min", c(rep(1,2*n),0), rbind(A1, A2),c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,1) [1] 1.077415 |
It looks like it’s working well…
Of course, we can adapt our previous code for quantiles
tau = .3 quantile(x,tau) 30% 0.6741586 |
The linear program is now\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbracewith a_i,b_i\geq 0 and y_i-\mu=a_i-b_i, \forall i=1,\cdots,n. The R code is now
A1 = cbind(diag(2*n),0) A2 = cbind(diag(n), -diag(n), 1) r = lp("min", c(rep(tau,n),rep(1-tau,n),0), rbind(A1, A2),c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,1) [1] 0.6741586 |
So far so good…
Consider the following dataset, with rents of flat, in a major German city, as function of the surface, the year of construction, etc.
base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE) |
The linear program for the quantile regression is now\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbracewith a_i,b_i\geq 0 and y_i-[\beta_0^\tau+\beta_1^\tau x_i]=a_i-b_i\forall i=1,\cdots,n. So use here
require(lpSolve) tau = .3 n=nrow(base) X = cbind( 1, base$area) y = base$rent_euro A1 = cbind(diag(2*n), 0,0) A2 = cbind(diag(n), -diag(n), X) r = lp("min", c(rep(tau,n), rep(1-tau,n),0,0), rbind(A1, A2), c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,2) [1] 148.946864 3.289674 |
Of course, we can use R function to fit that model
library(quantreg) rq(rent_euro~area, tau=tau, data=base) Coefficients: (Intercept) area 148.946864 3.289674 |
Here again, it seems to work quite well. We can use a different probability level, of course, and get a plot
plot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")), ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5) sf=0:250 yr=r$solution[2*n+1]+r$solution[2*n+2]*sf lines(sf,yr,lwd=2,col="blue") tau = .9 r = lp("min", c(rep(tau,n), rep(1-tau,n),0,0), rbind(A1, A2), c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,2) [1] 121.815505 7.865536 yr=r$solution[2*n+1]+r$solution[2*n+2]*sf lines(sf,yr,lwd=2,col="blue") |
Now that we understand how to run the optimization program with one covariate, why not try with two ? For instance, let us see if we can explain the rent of a flat as a (linear) function of the surface and the age of the building.
require(lpSolve) tau = .3 n=nrow(base) X = cbind( 1, base$area, base$yearc ) y = base$rent_euro A1 = cbind(diag(2*n), 0,0,0) A2 = cbind(diag(n), -diag(n), X) r = lp("min", c(rep(tau,n), rep(1-tau,n),0,0,0), rbind(A1, A2), c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,3) [1] 0.000000 3.257562 0.077501 |
Unfortunately, this time, it is not working well…
library(quantreg) rq(rent_euro~area+yearc, tau=tau, data=base) Coefficients: (Intercept) area yearc -5542.503252 3.978135 2.887234 |
Results are quite different. And actually, another technique can confirm the later (IRLS – Iteratively Reweighted Least Squares)
eps = residuals(lm(rent_euro~area+yearc, data=base)) for(s in 1:500){ reg = lm(rent_euro~area+yearc, data=base, weights=(tau*(eps>0)+(1-tau)*(eps<0))/abs(eps)) eps = residuals(reg) } reg$coefficients (Intercept) area yearc -5484.443043 3.955134 2.857943 |
I could not figure out what went wrong with the linear program. Not only coefficients are very different, but also predictions…
yr = r$solution[2*n+1]+r$solution[2*n+2]*base$area+r$solution[2*n+3]*base$yearc plot(predict(reg),yr) abline(a=0,b=1,lty=2,col="red") |
It’s now time to investigate….
Le schéma – ou diagramme – des relations entre les modes convergences ressemble à ce qui est présenté sur la droite. C’est grosso modo ce que j’ai dit en cours (sans le représenter explicitement).
Je ne vais pas démontrer ces implications car ce sont des théorèmes (qui peuvent donc se trouver dans n’importe quel textbook de probabilité).
En revanche, un étudiant m’avait demandé si les réciproques étaient vraies. Comme toujours en maths, si les réciproques étaient toujours vraies, on le mettrait. Donc ici les réciproques sont généralement fausses. On peut toutefois se demander quand ces réciproques peuvent être vraies…et quand elles sont fausses.
Commençons par le plus simple (à mon avis). Il existe un cas particulier où cette réciproque est vraie: le cas de variables aléatoires bornées. Autrement dit, si
et que l’on a
pour tout
, alors
pour tout
.
Bon, en fait la condition est un peu trop forte, et il suffit que la suite soit uniformément intégrable.
Maintenant, il existe aussi des cas où cette réciproque sera fausse (et on se doute qu’en construisant une suite non bornée, on devrait y arriver). Par exemple, supposons que pour tout ,
On peut montrer que , car en effet
Toutefois, certes la suite converge dans L1, mais vers une limite différente, étant donné que pour tout
. De plus,
et donc la suite ne converge pas dans L2.
Un exemple simple ici: si
converge en loi vers une constante, alors on aura convergence en probabilité (vers la même constante).
Maintenant la réciproque est relativement simple, l’idée étant que n’implique en rien
, au sens où
, pour tout
. Supposons maintenant que l’on ait à la fois
et
, où on va supposer que
est une loi symétrique. La symétrie signifie que
. Notons que l’on a alors forcément
. Supposons que la réciproque soit toujours vraie, alors on devrait avoir
, et donc on aurait que
, ce qui n’arrive que si
est la constante
.
Je ne vais pas m’étendre sur la construction d’exemples pour lesquels la réciproque est fausse (formellement, il suffit de changer d’espace de probabilité sous-jacent, je peux renvoyer vers les notes de cours de Philippe Briand (ici) et Hélène m’a suggéré également les exercices en ligne là). Par contre, il existe des cas importants où la réciproque est vraie.
Considérons une suite de variables indépendantes et identiquement distribuées . Posons alors
On peut alors montrer que
si et seulement si
(ça doit être du Borel-Cantelli).
Je ne vais pas non plus m’étendre davantage, mais on a un résultat qui ressemble à ce que je disais auparavant. En fait, il existe un théorème (dit de Skorohod) qui permet de dire que si on a une convergence en loi, alors on peut construite une suite de variables aléatoires qui
convergera presque sûrement. Formellement on a que si , alors il existe une suite
et une variable
, telles que
et
telles que
. Mais de manière générale, si la suite n’est pas construite comme il faut (
par
en quelques sortes) alors il est assez rare d’avoir la réciproque.