Monthly Archives: October 2012

ACT6420, Computer, Regression

Why pictures are so important when modeling data?

1 Comment

(bis repetita) Consider the following regression summary,

Call:
lm(formula = y1 ~ x1)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0001     1.1247   2.667  0.02573 *
x1            0.5001     0.1179   4.241  0.00217 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared: 0.6665,	Adjusted R-squared: 0.6295
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217
obtained from
> data(anscombe)
> reg1=lm(y1~x1,data=anscombe)
Can we say something if we look (only) at that output ? The intercept is significatively non-null, as well as the slope, the  is large (66%). It looks like we do have a nice model here. And in a perfect world, we might hope that data are coming from this kind of dataset,
But it might be possible to have completely different kinds of patterns. Actually, four differents sets of data are coming from Anscombe (1973). And that all those datasets are somehow equivalent: the ‘s have the same mean, and the same variance
> apply(anscombe[,1:4],2,mean)
x1 x2 x3 x4
9  9  9  9
> apply(anscombe[,1:4],2,var)
x1 x2 x3 x4
11 11 11 11
and so are the ‘s
> apply(anscombe[,5:8],2,mean)
y1       y2       y3       y4
7.500909 7.500909 7.500000 7.500909
> apply(anscombe[,5:8],2,var)
y1       y2       y3       y4
4.127269 4.127629 4.122620 4.123249
Further, observe also that the correlation between the ‘s and the ‘s is the same
> cor(anscombe)[1:4,5:8]
y1         y2         y3         y4
x1  0.8164205  0.8162365  0.8162867 -0.3140467
x2  0.8164205  0.8162365  0.8162867 -0.3140467
x3  0.8164205  0.8162365  0.8162867 -0.3140467
x4 -0.5290927 -0.7184365 -0.3446610  0.8165214
> diag(cor(anscombe)[1:4,5:8])
[1] 0.8164205 0.8162365 0.8162867 0.8165214
which yields the same regression line (intercept and slope)
> cbind(coef(reg1),coef(reg2),coef(reg3),coef(reg4))
[,1]     [,2]      [,3]      [,4]
(Intercept) 3.0000909 3.000909 3.0024545 3.0017273
x1          0.5000909 0.500000 0.4997273 0.4999091
But there is more. Much more. For instance, we always have the standard deviation for residuals
> c(summary(reg1)$sigma,summary(reg2)$sigma,
+ summary(reg3)$sigma,summary(reg4)$sigma)
[1] 1.236603 1.237214 1.236311 1.235695
Thus, all regressions here have the same R2
> c(summary(reg1)$r.squared,summary(reg2)$r.squared,
+ summary(reg3)$r.squared,summary(reg4)$r.squared)
[1] 0.6665425 0.6662420 0.6663240 0.6667073
Finally, Fisher’s F statistics is also (almost) the same.
+ c(summary(reg1)$fstatistic[1],summary(reg2)$fstatistic[1],
+ summary(reg3)$fstatistic[1],summary(reg4)$fstatistic[1])
value    value    value    value
17.98994 17.96565 17.97228 18.00329
Thus, with the following datasets, we have the same prediction (and the same confidence intervals). Consider for instance the second dataset (the first one being mentioned above),
> reg2=lm(y2~x2,data=anscombe)
The output is here exactly the same as the one we had above
> summary(reg2)

Call:
lm(formula = y2 ~ x2, data = anscombe)

Residuals:
Min      1Q  Median      3Q     Max
-1.9009 -0.7609  0.1291  0.9491  1.2691

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    3.001      1.125   2.667  0.02576 *
x2             0.500      0.118   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared: 0.6662,	Adjusted R-squared: 0.6292
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179
Here, the perfect model is the one obtained with a quadratic regression.
> reg2b=lm(y2~x2+I(x2^2),data=anscombe)
> summary(reg2b)

Call:
lm(formula = y2 ~ x2 + I(x2^2), data = anscombe)

Residuals:
Min         1Q     Median         3Q        Max
-0.0013287 -0.0011888 -0.0006294  0.0008741  0.0023776

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.9957343  0.0043299   -1385   <2e-16 ***
x2           2.7808392  0.0010401    2674   <2e-16 ***
I(x2^2)     -0.1267133  0.0000571   -2219   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001672 on 8 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 7.378e+06 on 2 and 8 DF,  p-value: < 2.2e-16
Consider now the third one
> reg3=lm(y3~x3,data=anscombe)
i.e.
> summary(reg3)

Call:
lm(formula = y3 ~ x3, data = anscombe)

Residuals:
Min      1Q  Median      3Q     Max
-1.1586 -0.6146 -0.2303  0.1540  3.2411

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0025     1.1245   2.670  0.02562 *
x3            0.4997     0.1179   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared: 0.6663,	Adjusted R-squared: 0.6292
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176
This time, the linear model could have been perfect. The problem is one outlier. If we remove it, we have
> reg3b=lm(y3~x3,data=anscombe[-3,])
> summary(reg3b)

Call:
lm(formula = y3 ~ x3, data = anscombe[-3, ])

Residuals:
Min         1Q     Median         3Q        Max
-0.0041558 -0.0022240  0.0000649  0.0018182  0.0050649

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.0056494  0.0029242    1370   <2e-16 ***
x3          0.3453896  0.0003206    1077   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.003082 on 8 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 1.161e+06 on 1 and 8 DF,  p-value: < 2.2e-16
Finally consider
> reg4=lm(y4~x4,data=anscombe)
This time, there is an other kind of outlier, in ‘s, but again, the regression is exactly the same,
> summary(reg4)

Call:
lm(formula = y4 ~ x4, data = anscombe)

Residuals:
Min     1Q Median     3Q    Max
-1.751 -0.831  0.000  0.809  1.839

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0017     1.1239   2.671  0.02559 *
x4            0.4999     0.1178   4.243  0.00216 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared: 0.6667,	Adjusted R-squared: 0.6297
F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165
The graph is here
So clearly, looking at the summary of a regression does not tell us anything… This is why we do spend some time on diagnostic, looking at graphs with the errors (the graphs above could be obtained only with one explanatory variable, while errors can be studied in any dimension): everything can be seen on thise graphs. E.g. for the first dataset,
or the second one
the third one
or the fourth one,
Computer, Maps, Maths, Miscellaneous

Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 2)

Leave a comment

Continuons la discussion de ce matin, avec l’exemple des capitales américaines. On peut trouver une telle liste liste en ligne sur wikipedia,

> capital=read.table(
+ "http://freakonometrics.free.fr/US-capital.csv",
+ header=TRUE,sep=";")
> capital=capital[-c(2,11),]

On va enlever Hawaï et l’Alaska, qui forcent à faire de trop grands détours. Dans TSP, sous R, il y a des listes de villes américaines, avec les coordonnées. En bricolant un peu, on va récupérer les coordonnées de la plupart des capitales,

> library(maps)
> library(maptools)
> library(sp)
> library(TSP)
> Lcap=paste(as.character(capital[,1]),
+ as.character(capital[,2]),sep=", ")
> M=as.matrix(USCA312)
> L=labels(M)[[1]]
> Lcap=Lcap[which(Lcap%in%L)]
> k=which(L%in%Lcap)
> M=M[k,]
> M=M[,k]
> listeUSA=TSP(M,L[k])

Et cette fois on est prêt,

> tour=solve_TSP(listeUSA, method = "nn")

on a lancé l’algorithme, et on peut faire un dessin,

> COORD=data.frame(USCA312_coords[k,])
> COORD=COORD[as.numeric(tour),] 
> map('state',fill=TRUE,col="yellow")
> lines(COORD$coords.x1,COORD$coords.x2,
+ lwd=3,col="blue")
> points(COORD$coords.x1,COORD$coords.x2,
+ pch = 19, cex = 1.2, col = "red")

pas mal, n’est pas ? Bon, mais comme d’habitude, rien ne garantie que l’on ait trouvé le trajet optimal… En fait, le soucis est que si on relance l’algorithme, on retombe sur un autre trajet,

voire encore un autre si on lancer une nouvelle fois l’algorithme,

C’est pénible, n’est ce pas ? En fait, si on lance 1,000 fois l’algorithme on obtient la distribution (des distances optimales) suivante

autrement dit, on a encore beaucoup de variabilité, avec une distance optimale qui peut varier d’environ 10% entre deux boucles d’algorithmes. Cela dit, cela reste faible par rapport à un trajet au hasard,

Et d’ailleurs, si on compare la distance obtenue en tirant les villes au hasard, on arrive à un taux de l’ordre de 27% (qui peut être comparé à ce que nous avions obtenu en tirant au hasard dans le carré unité)

>  HAZARD=OPTIMAL=RATIO=rep(NA,500)
>  for(s in 1:500){
+   HAZARD[s]=tsp.longueur(M,sample(1:nrow(M)))
+   tour=solve_TSP(listeUSA, method = "nn")
+   OPTIMAL[s]=tsp.longueur(M,as.numeric(tour))
+   RATIO[s]=HAZARD[s]/
+            OPTIMAL[s]
+ }

Bref, compte tenu de la (grande) dispersion des résultats optimaux obtenu, on pourrait dire que cet algorithme est assez mauvais. Même s’il est très rapide. Un autre avantage est que l’on peut spécifier la ville de départ, par exemple. Ici, si on souhaite partir de Floride, on peut obtenir le trajet suivant

> tour=solve_TSP(listeUSA, method = "nn",
+ control=list(start=41))

Passionnant n’est-ce pas ? Et encore, ça ne va pas m’aider à trouver le parcours optimal pour les enfants, car pour les enfants, il faut spécifier le point de départ (voire le point d’arrivée), il faut en plus rester sur la route (voire traverser aux passages pour piétons). On peut aussi vouloir faire une pause pipi au milieu de la promenade…. Bref, on a pas mal de contraintes ! Mais tout est expliqué dans le livre fabuleux de William Cook, In Pursuit of the Traveling Salesman, qui reprend toute l’histoire des algorithmes les plus fins pour résoudre ce joli problème de maths…

Computer, Maps, Maths, Miscellaneous

Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 1)

2 Comments

Bon, ben cette fois ça y est, c’est l’Halloween… et comme tous les ans, mon rôle devrait se limiter à creuser les citrouilles, décorer la maison… et attendre à la porte de la maison que les enfants des voisins viennent me dévaliser ma réserve de bonbons ! Je pourrais à la rigueur revendiquer un petit quelque chose dans le déguisement de mon fils (on a passé pas mal de temps à feuilleter les premiers tomes de Walking Deads, histoire de mieux saisir l’essence des personnages de morts-vivants, en plus de la Zombie Walk la fin de semaine passée, en centre-ville).

Mais cette année, pas de billet sur les zombies comme l’an passé, mais un vrai problème pratique: comment optimiser son trajet, quand on sait où sont les bonnes maisons, où on a de fortes chances d’avoir plein de bonbons (en particulier chez les parents des amis) ? C’est assez proche de l’idée de passer par un maximum de maisons, et ce à une heure raisonnable (histoire qu’il reste des bonbons à récupérer).

Pour tous ceux qui ont fait un peu de programmation dynamique, c’est un problème classique, connu sous le nom du “voyageur de commerce“. Ce problème est né en juillet 1954 avec “the shortest route for a traveling salesman” paru dans Newsweek ci-dessous (que l’on retrouve en détails dans le livre 50 Years of Integer Programming)

Enfin, disons qu’il a été popularisé à cette époque, car Euler évoquait déjà un problème similaire quand il travaillait sur le problème des sept ponts de Königsberg, dans “Solutio problematis ad geometriam situs pertinentis“,

(même si ce problème consiste à étudier l’existence d’un trajet sous des contraintes de passages par des nœuds, et pas vraiment de l’optimisation). On retrouve aussi un exemple en 1832, en Allemagne, dans un livre intitulé Der Handlungsreisende – wie er sein soll und was er zu thun hat, um Auftraege zu erhalten und eines gluecklichen Erfolgs in seinen Geschaeften gewiss zu sein – Von einem alten Commis-Voyageur où un trajet passant par 45 villes en Allemagne et en Suisse devait être optimisé

(le schéma ci-dessus montre celle décrite en 1832, pour un total de 1,285km, mais il est possible de faire une boucle de 1,248km).

En fait le problème a explosé en 1962 (à cette époque, il semble que les gens pouvaient encore penser que faire des maths pouvait être le fun), lors que Procter & Gamble ont offert un prix de 10,000$ à celui qui trouverait le chemin le plus court pour que Toody et Muldoom, les conducteurs de la voiture 54 dans une série populaire à l’époque, passent par 33 villes, aux États-Unis.

Cela dit, Carl Menger avait noté dès 1930 a noté que ce problème devrait etre difficile à résoudre. Facile à écrire, mais difficile à résoudre explicitement. Dans son livre, William Cooke rappelle qu’avec 33 villes, il y aurait

131,565,418,466,846,765,083,609,006,080,000,000

chemins possibles, dont il faudrait calculer la longueur. Et qu’avec l’ordinateur le plus puissant en 2009 (effectuant 1.5 million de milliards d’opérations à la seconde) il faudrait 28,000 milliards d’année pour mener à bien les calculs. Même si on essayait de tenir compte des progrès de l’informatique dans les mille prochaines années, ça risque de prendre du temps !..

En 1967, Jack Edmonds notait “I conjecture that there is no good algorithm for the traveling salesman problem” où “good” est une manière élégante pour dire un temps polynomial en , le nombre de villes. Ce qui s’appelle les problèmes de la classe . Il semblerait que le problème du voyageur de commerce soit un problème  et pas de la classe (cf le chapitre du livre de Sanjeev Arora et Boaz Barak, ou le chapitre 11 des notes de cours de Olivier Bournez, que j’ai pu découvrir grâce à @dmonniaux, ou formellement de la classe  et pas , i.e.  complet), voir aussi la discussion sur wikipedia),

Et depuis, l’institut Clay a proposé $1,000,000 à celui qui prouverait que les problèmes   ne sont pas de type  (ou qui prouveraient qu’ils le sont… laconjecture reste ouverte, semble-t-il). Bref, les enjeux financiers derrière ont considérablement augmenté (on pourra consulter Sipser (2007) pour la petite histoire)

Bon, et si l’histoire est intéressante, ce n’est pas une raison pour la lire en se croisant les bras… essayons de programmer un algorithme (simple) afin d’approcher une solution, car je ne pourrais me contenter de dire à mes enfants “c’est un problème  complet, papa a capitulé“.

Commençons par tirer 5 points au hasard dans le carré unité,

>  n = 5
>  (x = matrix(runif(2*n),nr=n))
[,1]       [,2]
[1,] 0.4142621 0.35600725
[2,] 0.9477231 0.66938653
[3,] 0.2518656 0.19900873
[4,] 0.5314231 0.22954670
[5,] 0.5156522 0.09066741
et commençons par calculer la distance entre les points (i.e. une matrice, ou ici une matrice triangulaire inférieure)
>  (d = dist(x))
1         2         3         4
2 0.6186980
3 0.2258786 0.8399243
4 0.1723919 0.6056111 0.2812205
5 0.2840514 0.7222195 0.2851688 0.1397719
On va mettre ce triangle sous forme de matrice (pleine), ça sera plus simple à manipuler par la suite
>  (d = as.matrix(d))
1         2         3         4         5
1 0.0000000 0.6186980 0.2258786 0.1723919 0.2840514
2 0.6186980 0.0000000 0.8399243 0.6056111 0.7222195
3 0.2258786 0.8399243 0.0000000 0.2812205 0.2851688
4 0.1723919 0.6056111 0.2812205 0.0000000 0.1397719
5 0.2840514 0.7222195 0.2851688 0.1397719 0.0000000
Le principe, pour initialiser l’algorithme est de tirer au hasard un trajet, c’est à dire une permutation de nos points,
>  (o = sample(1:n))
[1] 2 4 1 5 3
Il faut alors calculer le 4 distances entre ces 5 villes (on ne prévoit pas de retour à la ville de départ, pour l’instant). Le plus simple est de jouer avec cette matrice de la réarranger,
>  d[o[1:(n-1)],][,o[2:n]]
4         1         5         3
2 0.6056111 0.6186980 0.7222195 0.8399243
4 0.0000000 0.1723919 0.1397719 0.2812205
1 0.1723919 0.0000000 0.2840514 0.2258786
5 0.1397719 0.2840514 0.0000000 0.2851688
et de prendre les 4 valeurs sur la diagonale
>  diag(d[o[1:(n-1)],][,o[2:n]])
[1] 0.6056111 0.1723919 0.2840514 0.2851688
Ce sont les 4 distances que l’on cherche, entre les 5 villes. La distance totale parcourue est alors
>  sum(diag(d[o[1:(n-1)],][,o[2:n]]))
[1] 1.347223
L’idée est ensuite de permuter deux des villes
>  (i=sample(1:n,2))
[1] 3 1
(ici la première et la troisième de la liste, mais on pourrait décider de fixer la première si on doit partir de là) et de voir si la nouvelle distance est plus courte, ou pas,
>  os=o
>  os[i]=o[rev(i)]
>  sum(diag(d[os[1:(n-1)],][,os[2:n]]))
[1] 1.785391
La distance étant plus longue, on ne va pas permuter. On pourrait tenter 1000 fois de permuter deux villes. Si on ne parvient pas à réduire la distance totale, c’est qu’on a atteint une valeur intéressante (à défaut de pouvoir prouver qu’elle est optimale). Et si on peut améliorer la distance, on permute. Le plus simple est donc de faire des petites fonctions. La première reprend le calcul de la distance totale, sur la diagonale de la matrice permutée, à matrice de distance donnée et à permutation donnée,
> tsp.longueur=function(matrice,ordres) {
+ n=length(ordres)
+ sum(diag(matrice[ordres[1:(n-1)],][,ordres[2:n]])) }
 La seconde est juste une boucle,
> tsp.optimal=function(d,N=1000) {
+ d=as.matrix(d)
+ n=ncol(d)
+ o=sample(1:n)
+ v=tsp.longueur(d,o)
+ k=0
+ while(k < N) {
+   i<-sample(1:n,2)
+   os=o
+   os[i]=o[rev(i)]
+   w=tsp.longueur(d,os)
+   if(w < v) {
+     v=w
+     o=os
+     k=0} else {k=k+1}}
+ list(ordre=o,longueur=v)}
On peut se lancer maintenant, avec par exemple 30 points simulés, sur le carré unité,
>  set.seed(1)
>  n=30
>  x=matrix(runif(2*n),nr=n)
>  o=sample(1:n)
>  plot(x,xlim=0:1,ylim=0:1,xlab="",ylab="")
>  lines(x[o,],col='blue')
>  r=tsp.optimal(dist(x))
>  os=r$ordre
>  lines(x[os,],col='blue',lwd=1.5)
avec la première trajectoire, à gauche, qui reste en trait fin sur le dessin de droite,

Pas mal, non ? Bon, en fait le soucis est que si on joue plusieurs fois avec notre fonction, on obtient toujours un graphique différent,

Autrement dit, je trouve (certes) des chemins pas trop long, mais je suis loin de trouver le chemin optimal. Et histoire de clotûrer le débat, notons qu’on peut faire un peu de simulation, afin de comparer la distance obtenue en faisant le trajet au hasard, et un trajet “optimal“,

> RATIO=HASARD=OPTIMAL=matrix(NA,500,10)
> for(m in 1:10){
+ n=5*m
+ for(s in 1:500){
+  x=matrix(runif(2*n),nr=n)
+  r=tsp.opt(dist(x))
+  HAZARD[s,m]=tsp.longueur(as.matrix(dist(x)),1:n)
+  OPTIMAL[s,m]=r$longueur
+  RATIO[s,m]=HAZARD[s,m]/
+             OPTIMAL[s,m]
+ }
+}

Si on compare la longueur moyenne du trajet (sur 500 scénarios) avec un tour au hasard (en rouge) et ou tour optimisé (en bleu), on a

ce qui donne un ratio moyen (ou une moyenne des ratios, les deux sont représentés)

Autrement dit, avec une dizaine de points éparpillés (au hasard) dans le carré unité, un tour au hasard ne sera que deux fois plus long… étonnant non ?

Mais on doit pouvoir aller plus loin, parce que sous R, il y a des fonctions (et des packages) pour faire des algorithmes plus compliqués… à suivre avant la tournée des bonbons de demain soir !

ACT6420, Computer, Regression

Le passage au log dans les modèles linéaires

1 Comment

Un billet rapide pour compléter et illustrer le passage au log dans un modèle linéaire (que l’on abordera cette semaine en cours). Le point de départ est le modèle linéaire, où on suppose que, conditionnellement à ,  suit une loi normale. Pour rappel, si on a une loi normale, , alors  et . Les intervalles de confiance à 90% et 95% sont symétriques par rapport à la moyenne (qui est aussi la médiane, soit dit en passant),

Dans un modèle Gaussien avec homoscédasiticité,  i.e.  alors que . On a alors les bandes de confiance suivantes, pour un modèle de régression linéaire,

Bon, maintenant, que se passe-t-il si on prend l’exponentiel ? Pour la loi normale, rappelons que l’on obtient une loi lognormale, i.e. , les deux paramètres étant liés à la loi normale sous jacente, car désormais

alors que

Graphiquement, on a la loi suivante, avec les intervalles de confiance à 90% et 95% représentés ci-dessous. Le point noir est   alors que le point bleu est l’espérance de la loi lognormale.

On notera que le quantile de la loi log-normale est l’exponentiel du quantile de la loi normale. En effet, si https://latex.codecogs.com/gif.latex?\mathbb{P}(Y\leq%20q)=\alpha alors https://latex.codecogs.com/gif.latex?\mathbb{P}(\exp(Y)\leq%20\exp(q))=\alpha. En particulier,  n’est pas la moyenne de , mais la médiane (puisque  était la médiane de ).

Mais il n’est pas rare de voir utilisé un intervalle de confiance de la forme

qui est la forme classique de l’intervalle de confiance Gaussien (symétrique autour de la moyenne). Ici, on aurait les niveaux suivants

Notons qu’il n’y a aucune raison ici d’avoir une probabilité https://latex.codecogs.com/gif.latex?1-\alpha d’être dans l’intervalle de confiance obtenu avec les quantiles https://latex.codecogs.com/gif.latex?q_{1-\alpha/2} de la loi normale.

Maintenant, si on prend l’exponentiel d’un modèle linéaire (i.e. le logarithme de la variable d’intérêt est modélisé par un modèle linéaire) on a

avec une variance (conditionnelle) qui dépend de la variable explicative

Là encore, le plus naturel est d’utiliser comme bornes de l’intervalle de confiance des quantiles associés à la loi lognormale,

mais il n’est pas rare de voir utilisé des intervalles de type Gaussiens,

On perd là encore en interprétation car les bornes n’ont plus rien à voir avec les quantiles.

ACT2121, Course

ACT2121 premier examen

4 Comments

Cet après-midi avait lieu le premier examen d’actuariat 1, préparatoire pour l’examen P de la SOA. Comme pour l’examen P, il y avait des questions à choix multiples, en un temps (très limité). Le sujet est en ligne, ainsi que des éléments de corrections (bientôt) ici. Si vous n’êtes pas d’accord avec mes corrections, merci de me le dire rapidement, si possible avant que je n’attribue les notes…

ACT6420, Computer, Regression, Spatial

Une application du test de Fisher, partie 2

Leave a comment

Poursuivons un peu la discussion d’hier. En fait, si on poursuit, on arrive à regrouper des régions ensemble…. mais sans réelle cohérence car des régions proches au niveau électorale peuvent être éloignées géographiquement.

Si on veut que les régions regroupées soient proches, il faut une base avec des informations spatiales, comme savoir si deux départements se touchent, ou pas. Pour cela, on utilise la base suivante, qui contient

  • le numéro du département i
  • le numéro du département j
  • une variable indicatrice indiquant si les département se touche
  • la distance (en km) entre les centroïdes des départements
> france=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/departements-france.csv",
+ header=TRUE,sep=";")

On va nettoyer un peu, en enlevant les départements d’outre-mer,

> france0=france; i200=c(201,202)
> france0=france0[-c(which(france0$depi%in%i200),
+ which(france0$depj%in%i200)),]
> france0$ij=100*france0$depi+france0$depj

La base ressemble à ça,

> head(france0)
depi depj Cij Dij   ij
1   77   75   0  51 7775
2   78   75   0  39 7875
3   91   75   0  40 9175
4   92   75   1  10 9275
5   93   75   1  12 9375
6   94   75   1  13 9475

Plusieurs stratégies sont possible: regarde des départements voisins et voir si on peut les regrouper. Ou regarder des départements proches électoralement, et voir s’il sont proches géographiquement. Par exemple, si on ordonne les coefficients obtenus lors de la régression,

> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)
> COEF=coefficients(summary(reg))
> dbcoef=data.frame(COEF[,1:3],Dpt=rownames(COEF))
> rownames(dbcoef)=rank(COEF[,1])
> COEFtrie=dbcoef[as.character(1:94),]
> COEFtrie[1:6,]
Estimate Std..Error   t.value    Dpt
1 35.69903  0.3246129 109.97414 DptD06
2 36.46840  0.3244806 112.39010 DptD67
3 36.67594  0.3896297  94.13025 DptD68
4 37.39283  0.3202112 116.77552 DptD83
5 39.92023  0.3921762 101.79156 DptD74
6 42.61811  0.6182921  68.92877 DptD10

On peut se demander si les quatre premiers départements pourraient être regroupés ensemble. Une des modalités sert de modalité de référence, et on regarde les trois autres par rapport à cette dernière

> sbelection$Dpt=relevel(sbelection$Dpt,"D67")
> reg0=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> COEF0=coefficients(summary(reg0))
> COEF0[c(1,7,67,82),]
Estimate Std. Error     t value   Pr(>|t|)
(Intercept) 36.4684029  0.3244806 112.3900971 0.00000000
DptD06      -0.7693736  0.4589784  -1.6762740 0.09369320
DptD68       0.2075400  0.5070493   0.4093093 0.68231515
DptD83       0.9244243  0.4558759   2.0277980 0.04258819

On fait ensuite une analyse de la variance,

> sbelection$DptGroupe=sbelection$Dpt
> sblction$DptGroupe[sblction$Dpt%in%c("D06","D68","D83")]="D67"
> model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
> model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F   Pr(>F)
1  36114 2281072960
2  36111 2280191044  3    881917 4.6556 0.002954 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

On rejette ici l’hypothèse que les quatre départements se comportent de la même manière. Mais si on tente de regrouper les trois premières, cette fois, on accepte l’hypothèse

> sbelection$DptGroupe=sbelection$Dpt
> sblction$DptGroupe[sblction$Dpt%in%c("D06","D68")]="D67"
> model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
> model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36113 2280476723
2  36111 2280191044  2    285679 2.2621 0.1041

Sauf que sur les trois départements, seulement deux se touchent,

> france0[france0$ij%in%c(6768,667,668),]
depi depj Cij Dij   ij
3510    6   67   0 533  667
3552   67   68   1  95 6768
3605    6   68   0 439  668

On peut essayer de parcourir les sorties de la régression, et regarder si on a des voisins. Pour chaque département, on regarde si le voisin est un voisin géographique,

> departement=as.numeric(substr(COEFtrie$Dpt,5,6))
> deptij=departement[1:(nrow(COEFtrie)-1)]*100+
> departement[(2):(nrow(COEFtrie))]
> voisins=france0[france0$ij%in%deptij,]
> voisinage=tapply(voisins$Cij,factor(voisins$ij),max)
> COEFtrie$voisin=NA
> v=voisinage[as.character(deptij)]
> COEFtrie$voisin[1:(nrow(COEFtrie)-1)]=
> v[is.na(v)==FALSE]
> COEFtrie[1:10,]
Estimate Std..Error   t.value    Dpt voisin
1  35.69903  0.3246129 109.97414 DptD06      0
2  36.46840  0.3244806 112.39010 DptD67      1
3  36.67594  0.3896297  94.13025 DptD68      0
4  37.39283  0.3202112 116.77552 DptD83      0
5  39.92023  0.3921762 101.79156 DptD74      0
6  42.61811  0.6182921  68.92877 DptD10      0
7  42.77423  0.4397101  97.27826 DptD01      0
8  43.53395  0.4448212  97.86843 DptD84      0
9  44.38759  0.3897553 113.88579 DptD85      0
10 44.65361  0.4577710  97.54575 DptD51      1

On note qu’il y a plusieurs voisins. On peut ensuite faire un test de Fisher. S’il est positif, on fusionne les classes. Pour la régression, mais aussi dans la base des distances et des voisinages,

>  if(COEFtrie$voisin[i]==1){
+  sbelection$Dpt[sbelection$Dpt==substr(
+ COEFtrie$Dpt[i+vsn],4,6)]=substr(COEFtrie$Dpt[i],4,6)
+  france0$depi[france0$depi==as.numeric(
+  substr(COEFtrie$Dpt[i+vsn],5,6))]=as.numeric(
+ substr(COEFtrie$Dpt[i],5,6))
+  france0$depi[france0$depj==as.numeric(
+ substr(COEFtrie$Dpt[i+vsn],5,6))]=as.numeric(
+ substr(COEFtrie$Dpt[i],5,6))

Dans la base des départements, il y aura plusieurs paires. Si A et B ont fusionné, il suffit que C touche A ou B pour que C soit voisin des régions fusionnées. C’est pour cela que l’on utilise

> voisinage=tapply(voisins$Cij,factor(voisins$ij),max)

pour savoir si les départements sont voisins. On peut alors répéter plusieurs fois l’opération.

Par exemple, on note que 28, 45 et 89 peuvent fusionner,

> sbelection$Dpt=relevel(sbelection$Dpt,"D45")
+ reg10=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
+ COEF10=coefficients(summary(reg10))
+ COEF10[c(1,28,88),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 45.9303788  0.4170114 110.1417762 0.0000000
DptD28       0.5643212  0.6594028   0.8558065 0.3921105
DptD89       0.9423212  0.7016026   1.3430982 0.1792486
+ sbelection$DptGroupe=sbelection$Dpt
+ sbelection$DptGroupe[sbelection$Dpt%in%
+ c("D28","D89")]="D45"
+ model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
+ model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
+ anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36113 2280314007
2  36111 2280191044  2    122963 0.9737 0.3777

Géographiquement, ce sont effectivement des voisins

Mais peu de départements voisins peuvent effectivement être fusionnés.

ACT6420, Computer, Regression, Spatial

Une application du test de Fisher, partie 1

Leave a comment

Un petit (?) billet pour revenir sur un point évoqué mardi dernier en cours, sur l’utilisation du test de Fisher (et des techniques d’analyse de variance). Je vais juste évoquer rapidement un vieux billet sur l’histoire du 5% utilisé lorsque l’on calcule la -value (ce qui complètera la réponse faite en classe mardi lorsque la question a été posée). Mardi dernier, j’avais évoqué rapidement le test d’une contrainte linéaire sur les paramètres. On se demande ici si les deux modèles suivants sont vraiment différents, entre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-2.png

et sa version contrainte

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-1.png

On cherche alors à tester l’hypothèse

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-3.png

contre l’hypothèse alternative

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-4.png

Pour ça, on peut utiliser ce qu’a proposé Ronald Fisher, à savoir d’utiliser la somme des carrés des résidus (qui correspondent – à un facteur multiplicatif près – à l’estimateur de la variance du bruit), et il a posé

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-9.png

Je parlais tout à l’heure de choix de modèle… on notera que l’on en est pas très loin car on peut réécrire cette expression en faisant apparaître le  des régressions, contraintes et noncontraintes.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-10.png

Fisher a montré que sous l’hypothèse H0,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-11.png

i.e. la loi de la statistique de Fisher est précisément une loi de Fisher, Pour aller plus loin sur l’analyse de la variance (ANOVA, analysis of variance), je renvois àFaraway (2005).

J’avais dit oralement qu’une application intéressant du test de nullité de plusieurs coefficient dans une régression était son utilisation pour des variables catégorielles, et pour faire du regroupement de modalité.

Considérons un cas d’école. On dispose des résultats des votes dans plusieurs villes, en France, lors des dernières élections présidentielles (ces bases sont en ligne sur http://www.data.gouv.fr/), opposant François Hollande (avec un H) et Nicolas Sarkozy (avec un S),

> election=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/
205cf26f15a43974fdb7d3b25397613d.csv",
+ sep=";",header=TRUE,dec=",")
> head(election)
Code.du.dpt Inscrits Votants VoixH X.VoixExpH VoixS X.VoixExpS
1           1      592     506   195      41.05   280      58.95
2           1      215     175    82      52.23    75      47.77
3           1     8208    6447  3009      50.41  2960      49.59
4           1     1152     949   340      38.16   551      61.84
5           1      105      83    35      44.30    44      55.70
6           1     1702    1461   630      46.46   726      53.54

On va nettoyer un peu la base

> sbelection=election[election$Code.du.département
+ %in%as.character(1:95),]

car il y a des départements d’outre mers et la Corse, qui a une codification particulière. Bref, on retient les 95 départements métropolitain, moins le 20ème qui est la Corse. On va commencer par créer une variable de département qui soit un facteur qualitatif. Pour l’instant, le département est un numéro, entre 1 et 95 mais ce n’est aucunement un nombre !

> sbelection$Dpt=(paste("D",
> sbelection$Code.du.département,sep=""))
> I=as.numeric(as.character(sbelection$Code.du.département))<10
> sbelection$Dpt[I]=(paste("D0",
+ sbelection$Code.du.département,sep=""))
> sbelection$Dpt=factor(sbelection$Dpt)

On peut alors faire une régression, où on peut essayer de modéliser le taux de vote pour François Hollande (par exemple) en fonction de variables explicatives, et en particulier, on veut voir s’il n’y a pas un effet spatial. Bref, on va régresser sur le département. Et d’ailleurs que sur le département, inutile de mettre d’autre variable pour illustrer.

> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)

Par contre, je fais une régression pondérée, car les villes (et villages) ne font pas du tout la même taille. C’est exactement la même idée que de faire une moyenne pondérée pour avoir le taux moyen dans une région (on rediscutera des régressions pondérées en fin de cours).

> COEF=coefficients(summary(reg))

La sortie ressemble à

> COEF
Estimate Std. Error   t value Pr(>|t|)
DptD01 42.77423  0.4397101  97.27826        0
DptD02 52.39423  0.4572016 114.59765        0
DptD03 56.89434  0.5468760 104.03518        0
DptD04 51.06544  0.7807748  65.40355        0
DptD05 50.70525  0.8403979  60.33481        0
DptD06 35.69903  0.3246129 109.97414        0
DptD07 53.44185  0.5564489  96.04089        0
DptD08 51.86319  0.6360901  81.53434        0
DptD09 64.67228  0.8013389  80.70528        0
DptD10 42.61811  0.6182921  68.92877        0

[…]

DptD89 46.87270  0.5642231  83.07476        0
DptD90 50.47561  0.9063735  55.68964        0
DptD91 53.41748  0.3172310 168.38670        0
DptD92 49.54305  0.2850597 173.79888        0
DptD93 65.29062  0.3338954 195.54217        0
DptD94 56.48855  0.3196643 176.71211        0
DptD95 53.89382  0.3344373 161.14778        0

On peut d’ailleurs faire une prédiction, par département (en prenant comme prédiction celle obtenue pour la première ville du département)

> chg=which(diff(as.numeric(as.character(
+ sbelection$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)

On peut ensuite essayer de faire un petite dessin pour visualiser cette prédiction,

> chg=which(diff(as.numeric(as.character(
+ sbelection$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)
> nomdpt=read.table(
+ "http://freakonometrics.blog.free.fr/
public/data/departements.txt",sep=";")
> nomdpt$nom=tolower(nomdpt$V2)
> dpt=as.character(c(paste("0",
+ 1:9,sep=""),c(10:19,21:95)))
> noms=as.character(nomdpt[nomdpt$V3%in%dpt,2])
> couleurliste=rev(rainbow(n=100))
> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur=couleurliste[trunc((PredDept-35)*3)]

C’est un peu long (et compliqué) car il faut identifier le bon département sur la carte R (qui est codée en lettres). D’ailleurs la fonction suivante ne marche pas

> match=match.map(francemap,dpt)
> francemap$names[which(is.na(match.map(
+ francemap,dpt[35]))==FALSE)]
[1] "Indre-et-Loire" "Indre"

car en l’occurrence, seule l’Indre aurait du être obtenue. On va faire la correspondance à la main,

> match=rep(NA,length(francemap$names))
> for(i in 1:length(francemap$names)){
+ if(sum(dpt==francemap$names[i])>0)
+ {match[i]=which(dpt==francemap$names[i])}
+ }

Une fois les départements identifiés, on colorie,

> color=couleur[match]
> map(database="france", fill=TRUE, col=color)

C’est joli, mais comme toujours c’est un peu compliqué de faire des cartes…

A titre d’information, les couleurs correspondent aux valeurs suivantes (avec en abscisse, le taux de vote obtenu par François Hollande)

Bref, on a beaucoup de modalité car il y a 95 département. Mais on devrait pouvoir regrouper. Si on ne sait pas trop par où commencer, prenons un département au hasard, disons le 35 (i.e. Rennes),

> k=34

On peut ensuite regarder les estimations des coefficients, ou mieux, les intervalles de confiances des différents estimateurs, et voir ceux qui chevauche celui de Rennes (on va prendre des intervalles de confiance à 70% histoire d’avoir des chances que les deux coefficients soient vraiment proches)

> alpha=.7
> CI.COEF=confint(reg,level=alpha)

On peut visualiser ces intervalles, avec une bande horizontale pour Rennes,

> plot(c(1:19,21:95),COEF[,1],ylim=c(35,
+ 60),pch=3,cex=.5)
> segments(c(1:19,21:95),CI.COEF[,1],
+ c(1:19,21:95),CI.COEF[,2])
> lk=which((CI.COEF[,1]<CI.COEF[k,2])&
+ (CI.COEF[,2]>CI.COEF[k,1]))
> segments(c(1:19,21:95)[lk],CI.COEF[lk,1],
+ c(1:19,21:95)[lk],CI.COEF[lk,2],col="red",lwd=2)

et en rouge, les intervalles qui la chevauchent

Ça y est, on a une liste potentielle de six départements à regrouper… On va prendre comme modalité de référence le département qui nous intéresse, et voir si les coefficients des autres départements sont (ou pas) significatifs.

> sbelection$Dpt=relevel(sbelection$Dpt,"D35")
> reg35=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)
> COEF35=coefficients(summary(reg35))
> COEF35[c(1,lk[lk<k]+1,lk[lk>k]),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 55.66226653  0.3275407 169.9399801 0.0000000
DptD11       0.59043789  0.6272106   0.9413710 0.3465210
DptD32       0.96845567  0.7855154   1.2328920 0.2176241
DptD36      -0.03007518  0.7414655  -0.0405618 0.9676455
DptD62       0.56238714  0.4247897   1.3239189 0.1855384
DptD75      -0.06226653  0.4098013  -0.1519432 0.8792326
DptD81      -0.11970183  0.6067040  -0.1972986 0.8435950

On a confirmation que ces six départements sont non significativement (individuellement) différents de celui de référence (et donc on accepte l’hypothèse de nullité du coefficient). Mais peut-on pour autant les regrouper? C’est là que l’on peut faire un test de Fisher. Pour cela on fait deux régressions. La première est celle avec un regroupement de départements en une modalité unique (les sept sont désormais ensemble)

> sbelection$DptGroupe=sbelection$Dpt
> sbelection$DptGroupe[sbelection$Dpt%
+ in%substr(rownames(COEF[lk,]),4,6)]=paste("D",k+1,sep="")
> model0=lm(X.VoixExpH~DptGroupe,
+ weights=Votants,data=sbelection)

et la seconde est celle que l’ont vient de faire (qui sera notre hypothèse alternative, i.e. il ne faut pas regrouper les départements ensemble car ils sont différents – ou au moins l’un est différent)

> model1=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)

On peut alors faire le test de Fisher (en continuant à corriger de l’effet taille de la population, et en faisant des sommes pondérées)

> SCR0=sum(sbelection$Votants*residuals(model0)^2)
> SCR1=sum(sbelection$Votants*residuals(model1)^2)
> dl0=model0$df.residual
> dl1=model1$df.residual
> (F=(SCR0-SCR1)/(dl0-dl1)/SCR1*dl1)
[1] 0.917397

On obtient la statistique de Fisher, et le plus simple est alors de calculer la -value

> 1-pf(F,dl0-dl1,dl1)
[1] 0.4809411

On a fait (sans le savoir) une analyse de la variance, qui peut se faire plus simplement (?) en utilisant la fonction suivante

> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36117 2280538612
2  36111 2280191044  6    347568 0.9174 0.4809

Et on accepte ici l’hypothèse que ces sept départements peuvent être regroupé en un seul.

Enfin, si statistiquement, il y a du sens à regrouper les modalités, on peut se demander si ce regroupement de régions géographique à du sens,

> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur="red"
> match=match.map(francemap,dpt)
> color=couleur[match]
> map(dataLa première  estbase="france", fill=TRUE, col=color)

Les régions que l’on souhaite associer à Rennes sont fort éloignées, et les regrouper n’a peut-être pas de sens… (à suivre donc).

ACT6420, Regression

Régression sur des variables orthogonales

Leave a comment

Hier, on a évoqué en cours le théorème de Frisch-Waugh (évoqué sur ce blog il y a quelques mois). Derrière, il y a l’estimation de modèles sur des sous-ensembles de valeurs. Formellement, on a un modèle de la forme

https://latex.codecogs.com/gif.latex?Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=\boldsymbol{X}_1\boldsymbol{\beta}_1+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon
Les équations normales sont ici (cf transparents du cours, ou équations écrites au tableau) https://latex.codecogs.com/gif.latex?\boldsymbol{X}%27\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\boldsymbol{X}%27Y, que l’on peut écrire sous la forme d’un système,

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\boldsymbol{X}_1%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_1%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_1%27Y%20\\%20\boldsymbol{X}_2%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
Le point intéressant est que si les variables dans le premier ensemble et les variables dans le second sont orthogonales, i.e. https://latex.codecogs.com/gif.latex?\boldsymbol{X}_1%27\boldsymbol{X}_2%20=\boldsymbol{0} et https://latex.codecogs.com/gif.latex?\boldsymbol{X}_2%27\boldsymbol{X}_1%20=\boldsymbol{0}, alors les termes diagonaux du système disparaissent, et on retrouve les équations normales des régressions simples, i.e.

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\boldsymbol{X}_1%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1=\boldsymbol{X}_1%27Y%20\\%20\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
On a ici une conséquence du théorème de Frisch-Waugh: si les variables sont orthogonales, les estimateurs dans la régression multiple coïncident avec ceux obtenus avec des régressions simples.

Afin de mieux clarifier ce point, écrivons ce qui se passerait avec seulement deux variables. On va associer la constante à la première, i.e.

https://latex.codecogs.com/gif.latex?Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=(1%20\%20\boldsymbol{X}_1)\begin{pmatrix}\boldsymbol{\beta}_0%20\\%20\boldsymbol{\beta}_1\end{pmatrix}+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon
Dans ce cas, les équations normales sont

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}(1%20\%20\boldsymbol{X}_1)\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0%20\\%20\widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}Y%20\\%20\boldsymbol{X}_2%27(1\%20\boldsymbol{X}_1)%20\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0%20\\%20\widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
Bref, il faut que https://latex.codecogs.com/gif.latex?\boldsymbol%20X_1 et https://latex.codecogs.com/gif.latex?\boldsymbol%20X_2 soient orthogonales pour que le terme sur la seconde diagonale disparaisse. C’est ce qu’on retrouve dans la plupart des livres d’économétrie, il suffit de chercher le théorème de Frisch-Waugh dans l’index (même si le théorème de Frisch-Waugh va beaucoup plus loin puisqu’il explique quels ajustement faire pour lier les estimateurs dans une régression simple, et dans le cas multiple).

Mais est-ce un résultat purement théorique ou peut-on l’invoquer dans un cas pratique ? Autrement dit, si deux variables semblent orthogonales, que se passent-il numériquement ? Si je simule des variables explicatives indépendantes, les estimateurs coïncident-ils (comme le suggère la théorie) ?

Considérons un modèle linéaire avec deux variables explicative que l’on supposera indépendantes, mais ayant toutes deux un impact sur la variable explicative. Le plus simple est de considérer un triple de variables telles que

https://latex.codecogs.com/gif.latex?\boldsymbol{Z}=\begin{pmatrix}%20Y%20\\%20X_1\\X_2\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}%20\mu_Y%20\\%20\mu_{X_1}\\\mu_{X_2}\end{pmatrix},\begin{pmatrix}%201%20&%20\star%20&%20\star%20\\%20\star%20&%201&0\\\star&0&1\end{pmatrix}\right)
où https://latex.codecogs.com/gif.latex?\star désigne des valeurs non nulles. Par exemple,

> set.seed(1)
> n=100
> S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3)
> S
[,1] [,2] [,3]
[1,]  1.0  0.4  0.3
[2,]  0.4  1.0  0.0
[3,]  0.3  0.0  1.0
> M=c(1,3,6)
> Z=rmnorm(n,M,varcov=S)
> Y=Z[,1]; X1=Z[,2]; X2=Z[,3]
> base=data.frame(Y,X1,X2)
> cor(base)
Y          X1          X2
Y  1.0000000  0.37796580  0.27859381
X1 0.3779658  1.00000000 -0.05060539
X2 0.2785938 -0.05060539  1.00000000

On a simulé ici deux variables indépendantes, mais toutes deux corrélées avec la variable d’intérêt. On notera que numériquement

https://latex.codecogs.com/gif.latex?\frac{1}{n}\sum_{i=1}^n%20X_{1,i}X_{2,i}%20\neq%20\frac{1}{n}\sum_{i=1}^n%20X_{1,i}\cdot\frac{1}{n}\sum_{i=1}^n%20X_{2,i}

(autrement dit, en moyenne, on simulera des données indépendantes, mais dans ce cas précis, les données ne sont pas rigoureusement orthogonales).

> sum(X1*X2)/n
[1] 18.20194
> sum(X1)*sum(X2)/n^2
[1] 18.25097
Mais un test rapide nous pousse à accepter l’hypothèse d’indépendance entre les deux variables (sans trop d’hésitations).

> cor.test(X1,X2)

Pearsons product-moment correlation

data:  X1 and X2
t = -0.5016, df = 98, p-value = 0.6171
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2445923  0.1472766
sample estimates:
cor
-0.05060539

Si on regarde maintenant la régression,

> R=lm(Y~X1+X2,data=base)
> summary(R)

Call:
lm(formula = Y ~ X1 + X2, data = base)

Residuals:
Min       1Q   Median       3Q      Max
-1.87885 -0.57656  0.07957  0.52719  1.77407

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.71180    0.28246  16.682  < 2e-16 ***
X1           0.37232    0.08441   4.411 2.66e-05 ***
X2           0.25975    0.07755   3.350  0.00115 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7954 on 97 degrees of freedom
Multiple R-squared: 0.2317,	Adjusted R-squared: 0.2159
F-statistic: 14.63 on 2 and 97 DF,  p-value: 2.803e-06

Quel est le lien avec les régressions simples, sur chacune des variables ? Si on régresse sur la première, on obtient

> R1=lm(Y~X1,data=base)
> summary(R1)

Call:
lm(formula = Y ~ X1, data = base)

Residuals:
Min       1Q   Median       3Q      Max
-1.78973 -0.54060  0.04036  0.44138  1.94389

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.03165    0.27934  18.013  < 2e-16 ***
X1           0.35802    0.08859   4.041 0.000106 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8358 on 98 degrees of freedom
Multiple R-squared: 0.1429,	Adjusted R-squared: 0.1341
F-statistic: 16.33 on 1 and 98 DF,  p-value: 0.0001058

et si on régresse sur la seconde

> R2=lm(Y~X2,data=base)
> summary(R2)

Call:
lm(formula = Y ~ X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-1.8652 -0.5507  0.1104  0.6109  2.1435

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.85053    0.12495  46.824  < 2e-16 ***
X2           0.24244    0.08443   2.872  0.00501 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.867 on 98 degrees of freedom
Multiple R-squared: 0.07761,	Adjusted R-squared: 0.0682
F-statistic: 8.246 on 1 and 98 DF,  p-value: 0.005005

Si les données ne sont pas rigoureusement orthogonales, alors https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 dans le modèle multiple est légèrement différent de https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 obtenu dans le modèle simple. Si on regarde la distribution de nos estimateurs (n’oublions pas qu’il y a une marge d’erreur), ils sont proches, ou au moins statistiquement proches (avec à gauche la distribution de https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 et à droite celle de https://latex.codecogs.com/gif.latex?\widehat{\beta}_2)

Si on génère 100 bases de données avec ce modèle (avec 500 observations à chaque fois) on obtient les estimateurs suivants (à gauche les https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 et à droite les https://latex.codecogs.com/gif.latex?\widehat{\beta}_2, avec en abscisse le modèle multiple et en ordonnées, le modèle simple).

(la valeur centrale étant la vraie valeur, utilisée pour les simulations). On pourrait se demander si les valeurs éloignées du barycentre (ou de cette première bissectrice) correspondent à des cas où l’hypothèse d’indépendance entre les deux variables est rejetée ? On peut mettre des couleurs différentes en fonction de la https://latex.codecogs.com/gif.latex?p-value: plus la https://latex.codecogs.com/gif.latex?p-value est faible (rejet de l’hypothèse d’indépendance) plus la couleur sera vive (et rouge), et plus la https://latex.codecogs.com/gif.latex?p-value est forte (acceptation de l’hypothèse d’indépendance) plus la couleur sera terne (et jaune, voire blanche),

sur 1,000 simulations, on obtient

Le code est ici

> set.seed(1)
> library(mnormt)
> S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3)
> BX1=BX2=matrix(NA,1000,2)
> P=rep(NA,1000)
> for(s in 1:1000){
+   n=500
+   Z=rmnorm(n,c(6,3,1),varcov=S)
+   Y=Z[,1]; X1=Z[,2]; X2=Z[,3]
+   base=data.frame(Y,X1,X2)
+   R=lm(Y~X1+X2,data=base)
+   BX1[s,1]=coefficients(R)[2]
+   BX2[s,1]=coefficients(R)[3]
+   BX1[s,2]=coefficients(lm(Y~X1,data=base))[2]
+   BX2[s,2]=coefficients(lm(Y~X2,data=base))[2]
+   P[s]=cor.test(X1,X2)$p.value
+ }
> COL=(heat.colors(100))
> plot(BX1[,1],BX1[,2],pch=19,cex=.6,col=COL[trunc(100*P)])
> abline(a=0,b=1,lty=2)

Autrement dit le fait que l’estimation dans le modèle multiple et les modèles simples diffèrent beaucoup n’a rien a voir avec la https://latex.codecogs.com/gif.latex?p-value d’un test d’indépendance entre les deux variables… étonnant non ? En revanche, si on rajoute une observation pour rendre le modèle parfaitement non-corrélé, par exemple,

> difference=sum(X1*X2)*(n+1)/sum(X1)-sum(X2)
> X1=c(X1,0)
> X2=c(X2,difference)
> n=n+1
> sum(X1*X2)/n
[1] 18.02172
> sum(X1)*sum(X2)/n^2
[1] 18.02172
On retrouve ici une corrélation parfaitement nulle,
> Y=c(Y,0)
>  base=data.frame(Y,X1,X2)
>  cor(base)
Y            X1            X2
Y  1.0000000  3.946091e-01  2.923412e-01
X1 0.3946091  1.000000e+00 -4.817519e-16
X2 0.2923412 -4.817519e-16  1.000000e+00
Si on regarde maintenant la régression multiple
>  R=lm(Y~X1+X2,data=base)
>  coefficients(R)
(Intercept)          X1          X2
-1.5028358   0.3589776   0.2531358
et les deux régressions simples,
>  R1=lm(Y~X1,data=base)
>  coefficients(R1)
(Intercept)          X1
0.02847362  0.35897764
>  R2=lm(Y~X2,data=base)
>  coefficients(R2)
(Intercept)          X2
-0.4334011   0.2531358
on obtient exactement les mêmes coefficient.
Moralité ? l’hypothèse d’orthogonalité entre les variables explicatives n’est pas à prendre à la légère: il faut une orthogonalité stricte pour que numériquement, la régression multiple soit équivalente aux régressions simples (comme le prédit la théorie). Mais de manière un peu surprenante, si on est très proche de l’indépendance, les valeurs numériques des estimateurs peut différer. Et ce, parfois de manière non négligeable, même si statistiquement, on valide l’hypothèse d’orthogonalité entre les variables explicatives.
2012, ACT2121

ACT2121, éléments de correction (2)

1 Comment

Un petit complément pour revenir sur un calcul fait hier, en classe, au tableau.

“Une police d’assurance couvre des pertes X et Y qui sont des variables aléatoires continues fonction dont on donne la densité conjointes. Trouver l’espérance du remboursement s’il est de X + Y avec un déductible de 1.”

Au tableau, je m’étais lancé dans des calculs un peu compliqué (notant que les variables étaient indépendantes, on pouvait utiliser la convolution pour trouver la loi de la somme facilement… mais on reviendra sur ce calcul plus tard). La solution la plus simple était la suivante: la loi de https://latex.codecogs.com/gif.latex?(X,Y) admet pour densité https://latex.codecogs.com/gif.latex?f_{X,Y}(x,y)=2xboldsymbol{1}_{[0,1]times[0,1]}(x,y). Le remboursement https://latex.codecogs.com/gif.latex?R est la variable aléatoire

https://latex.codecogs.com/gif.latex?R=h(X,Y)=(X+Y-1)boldsymbol{1}_{]1,+infty[}(X+Y)

On peut alors faire le calcul de
https://latex.codecogs.com/gif.latex?mathbb{E}(R)=int_{mathbb{R}^2}%20h(x,y)f(x,y)dxdy
i.e.
https://latex.codecogs.com/gif.latex?mathbb{E}(R)=int_0^1%20int_{1-x}^1%20[x+y-1]2x%20dy%20dx
soit
https://latex.codecogs.com/gif.latex?int_0^1%202xleft[(x-1)y+frac{y^2}{2}<br /> ight]_{1-x}^1%20dx=int_0^1%20x^3dx=frac{1}{4}
Voilà pour la version simple (et finalement assez générale). Pour la version plus compliquée
https://latex.codecogs.com/gif.latex?%20f_{X,Y}(x,y)=2xboldsymbol{1}_{[0,1]times[0,1]}(x,y)=underbrace{2x%20boldsymbol{1}_{[0,1]}(x)}underbrace{%20boldsymbol{1}_{[0,1]}(y)}
autrement dit, on a obtenu que les deux composantes étaient indépendantes, et que https://latex.codecogs.com/gif.latex?%20Y est uniformément distribuée sur https://latex.codecogs.com/gif.latex?%20[0,1]. Aussi, pour la loi de la somme
https://latex.codecogs.com/gif.latex?%20mathbb{P}(S%3Es)=int_{0}^1%20mathbb{P}(S%3Es|Y=y)f_Y(y)dy=int_{0}^1%20mathbb{P}(X%3Es-y|Y=y)dy
et par indépendance,

https://latex.codecogs.com/gif.latex?mathbb{P}(X%3Es-y|Y=y)=mathbb{P}(X%3Es-y)=int_{s-y}^1%20f_X(x)dx

i.e.
https://latex.codecogs.com/gif.latex?mathbb{P}(X%3Es-y|Y=y)=%20int_{s-y}^1%202xboldsymbol{1}_{[0,1]}(x)%20dx=int_{max{s-y,0}}^1%202xdx
soit

https://latex.codecogs.com/gif.latex?mathbb{P}(X%3Es-y|Y=y)=left[x^2<br /> ight]_{max{s-y,0}}^1=1-(max{s-y,0})^2

Aussi

https://latex.codecogs.com/gif.latex?mathbb{P}(S%3Es)=int_{0}^1%20big(%201-(max{s-y,0})^2%20big)dy

Première chose, on note que
https://latex.codecogs.com/gif.latex?mathbb{P}(S%3E1)=int_{0}^1%20big(%201-(max{1-y,0})^2%20big)dy=int_{0}^1%20big(%201-(1-y)^2%20big)dy
i.e.
https://latex.codecogs.com/gif.latex?mathbb{P}(S%3E1)=left[y+frac{(1-y)^3}{3}<br /> ight]_0^1=1+0-frac{1}{3}=frac{2}{3}
J’avais alors suggéré d’utiliser (en utilisant la propriété que pour une variable positive, l’espérance est l’intégrale de la fonction de survie) que
https://latex.codecogs.com/gif.latex?mathbb{E}(R)=int_0^1%20mathbb{P}(R%3Es)ds
Or cette expression n’est valide que pour une variable strictement positive ! Or ici, il y a de forte chance (1/3) que le remboursement soit nul. Le plus simple est alors de faire des calculs classiques d’espérance. Notons que
https://latex.codecogs.com/gif.latex?mathbb{P}(S%3E1+s)=left[y+frac{(1+s-y)^3}{3}%20<br /> ight]_0^1=1+frac{s^3}{3}-frac{(s+1)^3}{3}
donc un peu de calcul de dérivées permet d’écrire https://latex.codecogs.com/gif.latex?f_S(s)=s(2-s) pour https://latex.codecogs.com/gif.latex?sin[1,2]. On peut alors calculer l’espérance du remboursement
https://latex.codecogs.com/gif.latex?mathbb{E}(R)=int_1^2(s-1)f_S(s)ds=int_1^2(s-1)s(2-s)ds
soit
https://latex.codecogs.com/gif.latex?mathbb{E}(R)=int_1^2(s^3+3s^2-2s)ds=left[frac{s^4}{4}+s^3-s^2<br /> ight]_1^2=frac{1}{4}
Youpi, avec plein de calculs, on arrive au même résultat…

Faisons un peu de simulations pour vérifier,

> n=1000000 
> set.seed(1) 
> X=sqrt(runif(n)) 
> Y=runif(n)

On peut vérifier par exemple la densité de https://latex.codecogs.com/gif.latex?S sur https://latex.codecogs.com/gif.latex?[1,2],

> u=seq(1,2,by=.01) 
> hist(X+Y,probability=TRUE,col="light blue") 
> lines(u,u*(2-u),lwd=2,col="red")

Si on regarde maintenant la prime pure, on retrouve le résultat obtenu par deux calculs différents,

> R=(X+Y-1)*(X+Y>1) 
> mean(R) 
[1] 0.2497528
ACT2121

Actuariat 1, ACT2121, troisième cours

Leave a comment

Pour la troisième série d’exercices, les exercices sont en ligne ici.

Je rappelle que lundi 29 octobre aura lieu le premier examen (qui compte pour 33,33% de la note finale, comme annoncé dans le plan de cours, en ligne ici). Il y aura 30 questions, en français, du même type que celles vues toutes les semaines. La correction se fera suivant le principe de l’examen P de la SOA

“A candidate’s score will be based on the number of questions answered correctly. No credit will be given for omitted answers, and no credit will be lost for wrong answers; therefore, a candidate should answer all questions, even if the candidate needs to guess.“

Computer, Copulas, Geometry, Graphics, Maths

Fractals and Kronecker product

2 Comments

A few years ago, I went to listen to Roger Nelsen who was giving a talk about copulas with fractal support. Roger is amazing when he gives a talk (I am also a huge fan of his books, and articles), and I really wanted to play with that concept (that he did publish later on, with Gregory Fredricks and José Antonio Rodriguez-Lallena). I did mention that idea in a paper, writen with Alessandro Juri, just to mention some cases where deriving fixed point theorems is not that simple (since the limit may not exist).

The idea in the initial article was to start with something quite simple, a the so-called transformation matrix, e.g.

https://latex.codecogs.com/gif.latex?T=\frac{1}{8}\left(\begin{matrix}1&%200%20&%201%20\\%200%20&%204%20&%200%20\\%201%20&%200&1\end{matrix}\right)
Here, in all areas with mass, we spread it uniformly (say), i.e. the support of https://latex.codecogs.com/gif.latex?T(C^\perp) is the one below, i.e. https://latex.codecogs.com/gif.latex?1/8th of the mass is located in each corner, and https://latex.codecogs.com/gif.latex?1/2 is in the center. So if we spread the mass to have a copula (with uniform margin,)we have to consider squares on intervals https://latex.codecogs.com/gif.latex?[0,1/4], https://latex.codecogs.com/gif.latex?[1/4,3/4] and https://latex.codecogs.com/gif.latex?[3/4,1],

Then the idea, then, is to consider https://latex.codecogs.com/gif.latex?T^2=\otimes^2T, where  https://latex.codecogs.com/gif.latex?\otimes^2T is the tensor product (also called Kronecker product) of https://latex.codecogs.com/gif.latex?T with itself. Here, the support of https://latex.codecogs.com/gif.latex?T^2(C^\perp) is

Then, consider https://latex.codecogs.com/gif.latex?T^3=\otimes^3T, where https://latex.codecogs.com/gif.latex?\otimes^3T is the tensor product of https://latex.codecogs.com/gif.latex?T with itself, three times. And the support of https://latex.codecogs.com/gif.latex?T^3(C^\perp) is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if https://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}=(a_{i,j}), then

https://latex.codecogs.com/gif.latex?%20%20%20%20%20\mathbf{A}\otimes\mathbf{B}%20=%20\begin{pmatrix}%20a_{11}%20\mathbf{B}%20&%20\cdots%20&%20a_{1n}\mathbf{B}%20\\%20\vdots%20&%20\ddots%20&%20\vdots%20\\%20a_{m1}%20\mathbf{B}%20&%20\cdots%20&%20a_{mn}%20\mathbf{B}%20\end{pmatrix}

So, we need a transformation matrix: consider the following https://latex.codecogs.com/gif.latex?4\times4 matrix,

> k=4
> M=matrix(c(1,0,0,1,
+            0,1,1,0,
+            0,1,1,0,
+            1,0,0,1),k,k)
> M
[,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    1    0
[3,]    0    1    1    0
[4,]    1    0    0    1

Once we have it, we just consider the Kronecker product of this matrix with itself, which yields a https://latex.codecogs.com/gif.latex?4^2\times4^2 matrix,

> N=kronecker(M,M)
> N[,1:4]
[,1]  [,2] [,3] [,4]
[1,]     1    0    0    1
[2,]     0    1    1    0
[3,]     0    1    1    0
[4,]     1    0    0    1
[5,]     0    0    0    0
[6,]     0    0    0    0
[7,]     0    0    0    0
[8,]     0    0    0    0
[9,]     0    0    0    0
[10,]    0    0    0    0
[11,]    0    0    0    0
[12,]    0    0    0    0
[13,]    1    0    0    1
[14,]    0    1    1    0
[15,]    0    1    1    0
[16,]    1    0    0    1

And then, we continue,

> for(s in 1:3){N=kronecker(N,M)}

After only a couple of loops, we have a https://latex.codecogs.com/gif.latex?4^5\times4^5 matrix. And we can plot it simply to visualize the support,

> image(N,col=c("white","blue"))

As we zoom in, we can visualize this fractal property,

2012, ACT2121

ACT2121, éléments de correction

Leave a comment

La plupart des calculs pouvant se faire sans calculs (trop) complexes avec une calculatrice, je vais revenir sur un exercice (l’exercice 11 de la seconde série) pour proposer une correction.

“Le petit Nestor collectionne les cartes de joueurs de Baseball dans les paquets de gommes à mâcher. Il y a en tout 20 cartes différentes (réparties aléatoirement, une par paquet). Combien de paquets de gommes Nestor devrait-il s’attendre à avoir à acheter pour obtenir la collection complète ? ”

La bonne stratégie était d’ordonner les cartes par ordre d’apparition, pour la première fois, et de modéliser le nombre de paquets entre deux premières apparitions, comme indiqué sur le dessin ci-dessous


Si  désigne le nombre total de paquets à acheter, on note  le nombre de paquets à acheter, entre l’apparition de la ème nouvelle carte, et la ème (avec la convention que vaut 1, i.e. à l’achat du premier paquet, on a notre première carte). On notera que
https://latex.codecogs.com/gif.latex?M_i\sim\mathcal{G}((20-(i-1))/20)
et à partir de là, les calculs sont simples, puisque

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=mathbb{E}(M_1+M_2+cdots+M_{20})

soit (par linéarité de l’espérance)

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=mathbb{E}(M_1)+mathbb{E}(M_2)+cdots+mathbb{E}(M_{20})

i.e.

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=left(%201+%20frac{20}{19}%20+%20frac{20}{18}%20+%20frac{20}{17}%20+%20...%20+%20frac{20}{2}+frac{20}{1}%20<br /> ight)

ou encore

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=20left(%20frac{1}{20}+%20frac{1}{19}%20+%20frac{1}{18}%20+%20frac{1}{17}%20+%20...%20+%20frac{1}{2}+1%20<br /> ight)

La bonne réponse était alors (il faut sommer les vingt termes)

> sum(20/(1:20)) [1] 71.95479

Mais cette sommation de vingt termes n’est pas triviale, aussi, en cours, j’avais suggéré que

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)
> log(20)*20 [1] 59.91465

qui diffère du résultat numérique car il manquait la constante d’Euler

https://latex.codecogs.com/gif.latex?gamma%20=%20lim_{n%20<br /> ightarrow%20infty%20}%20left(%201+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20-%20ln(n)%20<br /> ight)

i.e.

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)+gamma

Numériquement, on obtient, en prenant pour https://latex.codecogs.com/gif.latex?gamma 0.57721 (cf Google)

> (log(20)+ 0.57721)*20 [1] 71.45885
2012, ACT2121

Actuariat 1, ACT2121, second cours

Leave a comment

Comme la dernière fois, je mets en ligne quelques exercices (une soixantaine cette fois). On devrait commencer à parler des lois usuelles de probabilités, pour l’instant des lois discrètes

  • la loi de Poisson

Si https://latex.codecogs.com/gif.latex?X suit une loi https://latex.codecogs.com/gif.latex?mathcal{P}(lambda),

https://latex.codecogs.com/gif.latex?%20%20%20%20%20p(k)%20=%20mathbb%20P(X%20=%20k)=%20e^{-lambda}frac{lambda%20^k}{k!},

Dans ce cas, https://latex.codecogs.com/gif.latex?mathbb{E}(X)=lambda et https://latex.codecogs.com/gif.latex?text{Var}(X)=lambda. La fonction génératrice des moments est alors

https://latex.codecogs.com/gif.latex?%20%20%20%20%20M_{X}(t)=%20mathbb{E}(e^{tX})=expleft(lambda%20(e^t-1)<br /> ight)

  • la loi binomiale

Si https://latex.codecogs.com/gif.latex?X suit une loi https://latex.codecogs.com/gif.latex?mathcal{B}(n,p),

https://latex.codecogs.com/gif.latex?%20%20%20%20%20p(k)%20=mathbb%20P(X%20=%20k)=%20{n%20choose%20k}%20,%20p^k%20q^{n-k}

Dans ce cas, https://latex.codecogs.com/gif.latex?mathbb%20E(X%20)=np et https://latex.codecogs.com/gif.latex?text{Var}(X%20)=np(1-p). La fonction génératrice des moments est alors

https://latex.codecogs.com/gif.latex?%20%20%20%20%20p(k)%20=%20mathbb%20P(X%20=%20k)=%20e^{-lambda}frac{lambda%20^k}{k!},

  • la loi binomiale négative

Enfin,, si https://latex.codecogs.com/gif.latex?X suit une loi https://latex.codecogs.com/gif.latex?%20BN(r,p),

 https://latex.codecogs.com/gif.latex?%20%20%20%20%20p(k)%20=%20mathbb%20P(X%20=%20k)=%20%20frac{Gamma(r+k)}{k!,Gamma(r)},p^r,q^k%20!

Dans ce cas,

https://latex.codecogs.com/gif.latex?%20mathbb%20E(X%20)=%20frac{r(1-p)}{p}!

et

https://latex.codecogs.com/gif.latex?%20text{Var}(X%20)=%20frac{r(1-p)}{p^2}!

La fonction génératrice des moments est alors

https://latex.codecogs.com/gif.latex?%20%20%20%20%20M_{X}(t)=%20mathbb{E}(e^{tX})=%20left(frac{p}{1-[1-p]%20e^t}<br /> ight)^r%20!

Je remettrais la semaine prochaine les lois continues les plus usuelles. Demain, soit on attaque cette nouvelle série, soit on revient sur des exercices de la semaine passée qui auraient posé problème…

Actuarial science, Computations, Computer

Compound Poisson and vectorized computations

3 Comments

Yesterday, I was asked how to write a code to generate a compound Poisson variables, i.e. a series of random variables  where  is a counting random variable (here Poisson disributed) and where the ‘s are i.i.d (and independent of ), with the convention  when . I came up with the following algorithm, but I was wondering if it was possible to get a better one…

>  rcpd=function(n,rN,rX){
+  N=rN(n)
+  X=rX(sum(N))
+  I=as.factor(rep(1:n,N))
+  S=tapply(X,I,sum)
+  V=as.numeric(S[as.character(1:n)])
+  V[is.na(V)]=0
+  return(V)}

Here, consider – to illustrate – the case where  and ,

>  rN.P=function(n) rpois(n,5)
>  rX.E=function(n) rexp(n,2)

We can generate a sample

>  S=rcpd(1000,rN=rN.P,rX=rX.E)

and check (using simulation) than 

> mean(S)
[1] 2.547033
> mean(rN.P(1000))*mean(rX.E(1000))
[1] 2.548309

and that 

> var(S)
[1] 2.60393
> mean(rN.P(1000))*var(rX.E(1000))+
+ mean(rX.E(1000))^2*var(rN.P(1000))
[1] 2.621376

If anyone might think of a faster algorithm, I’d be glad to hear about it…

ACT2040, ACT6420, Computer, Statistics, Tests

Test, valeur critique et p-value

14 Comments

Un petit complément suite au cours de mercredi dernier, pour insister sur l’importance de la p-value dans la lecture de la sortie d’un test.

  • Les erreurs dans un test statistique

Mais avant, rappelons qu’un test est une prise de décision: accepter ou rejeter une hypothèse. Et qu’on peut commettre une erreur. Ou pour être plus précis, on peut commettre deux types d’erreur,
• accepter l’hypothèse alors que cette dernière est fausse
• rejeter l’hypothèse alors que cette dernière était vraie
Pour reprendre une terminologie plus médicale, un test de grossesse peut dire à une femme qu’elle n’est pas enceinte, alors qu’elle l’est; ou dire qu’elle l’est, alors qu’elle ne l’est pas (voir tous les exemples dans les exercices de probabilités de l’examen P de la SOA, ou le cours ACT2121).
Formellement, on a deux probabilités,
• la probabilité d’accepter à tort notre hypothèse (on parlera d’erreur de second espèce), \beta
• la probabilité de rejeter à tort notre hypothèse (on parlera d’erreur de première espèce) \alpha
Dans un monde idéal on voudrait que les deux probabilités soient aussi petites que possibles… Mais c’est impossible, et le plus souvent, baisser une des probabilités se fait en augmentant l’autre. Les cas extrêmes étant
• avoir un test de grossesse qui déclare tout le monde enceinte: on ne rejette alors jamais à tort (on ne rejette jamais tout court en fait), mais on a un fort taux d’acceptation à tort,
• avoir un test de grossesse qui ne déclare personne enceinte: on n’accepte jamais à tort (car on n’accepte jamais) mais on a un fort taux de rejet à tort.
Bref, on a un arbitrage à faire entre deux types d’erreurs. Souvent, en pratique on va demander à contrôler l’erreur de première espèce (i.e. \alpha de l’ordre de 5%), et on chercher a un test qui, à \alpha donné, possède la plus faible erreur de première espèce. Voilà en gros pour la théorie: on se donne un seuil de significativité \alpha, qui correspond à la probabilité d’erreur de premier type. Et on va chercher à tester si une hypothèse H_0 est vraie, l’alternative étant une hypothèse H_1.

H_0 vraie H_1 vraie
accepter H_0 OK erreur
type 2
rejeter H_0 erreur
type 1
 OK
  • La “valeur critique”

La notion de valeur critique a été introduite dans Neyman & Pearson (1928). Cette valeur dépend de la forme de l’hypothèse alternative, en particulier savoir si le test est bilatéral, unilatéral à gauche, ou unilatéral à droite. Pour un test donné, la valeur critique peut-être vue comme la valeur limite a partir de laquelle on pourra rejeter H_0 avec un seuil de significativité donné.

  • La p-value

La p-value a été introduite dans Gibbons & Pratt (1975), meme si on peut retrouve l’idée beaucoup plus tôt, comme Pearson (1900), qui propose de calculer “the probability that the observed value of the chi-square statistic would be exceeded under the null hypothesis“. La p-value est la probabilité, sous H_0, d’obtenir une statistique aussi extrême (pour ne pas dire aussi grande) que la valeur observée sur l’échantillon. Aussi, pour un seuil de significativité \alpha donné, on compare p et \alpha, afin d’accepter, ou de rejeter H_0,
• si p\leq\alpha, on va rejeter l’hypothèse H_0 (en faveur de H_1)
• si p>\alpha, on va rejeter H_1 (en faveur de H_0).
On peut alors interpréter la p-value comme le plus petit seuil de significativité pour lequel l’hypothèse nulle est acceptée. Gibbons & Pratt (1975) reviennent longuement sur les interprétations, et surtout les mauvaises interprétations, de cette p-value.

  • Valeur critique versus p-value

Si on formalise un peu, on peut vouloir tester H_0:\theta=\theta_0 contre H_1:\theta>theta_0 (par exemple). De manière très générale, on dispose d’une statistique de test T qui a pour loi, sous H_0, F_{\theta_0}(\cdot) (que l’on supposera continue). Notons qu’on peut considérer une hypothèse alternative de la forme H_1:\theta\neq\theta_0, c’est juste plus pénible parce qu’il faut travailler sur \vert T\vert, et calculer des probabilités à gauche, ou à droite. Donc pour notre exemple, on va prendre un test unilatéral.
Dans l’approche classique (telle que présentée dans tous les cours de statistiques), on se donne un seul d’acceptation \alpha petit (disons 5%), et on cherche une valeur critique T_{1-alpha} telle que

Pour ceux qui se souviennent de leur cours de stats, cela peut faire penser à la puissance du test, définie par

\pi(\theta\vert \alpha)=\mathbb{P}(T\geq T_{1-\alpha}\vert \theta)=1-F_{\theta}(T_{1-\alpha})

Formellement, la p-value associée au test T est la variable aléatoire P définie par
P=1-F_{\theta_0}(T).
Donc effectivement, la p-value et la puissance sont liées, puisque

\mathbb{P}(P\leq \alpha\vert \theta)=\pi(\theta\vert \alpha)

autrement dit, la puissance peut-être vue comme la fonction de répartition de la p-value.

  • Intérêt computationnel de la p-value

D’un point de vue computationnel, la p-value est l’outil le plus important pour interpréter la sortie d’un test. Commençons par un test simple, comme une comparaison de moyennes. On cherche ici à tester H_0:\mu_X=\mu_Y contre H_1:\mu_X>\mu_Y pour des moyennes calculées sur deux groupes. Pour reprendre l’exemple abordé dans un précédant billet, on a les notes obtenues en ACT6420 par deux groupes différents. Et on veut savoir s’ils sont vraiment différents (ci-dessous le nombre de bonnes réponses, sur 40 questions, on travaillera ensuite sur la note sur 100)

image manquante

La statistique de test est ici

T = \frac{\overline{X} - \overline{Y}}{\displaystyle{ \sqrt{ {s_X^2 \over n_X} + {s_Y^2 \over n_Y} }}}

et sous H_0, T va suivre une loi de Student à \nu degrés de liberté, où \nu est donné par la relation de Welch–Satterthwaite (d’après Satterwaite (1946) et Welch (1947)),

\nu = {{\left( {s_X^2 \over n_X} + {s_Y^2 \over n_Y}\right)^2 } \over {{s_X^4 \over n_X^2 \cdot \left({n_X-1}\right)}+{s_Y^4 \over n_Y^2 \cdot \left({n_Y-1}\right)}}}

Numériquement, on a ici

> Xbar=mean(X)
> Ybar=mean(Y)
> Sx2=var(X)
> Sy2=var(Y)
> nX=length(X)
> nY=length(Y)
> (T=(Xbar-Ybar)/sqrt(Sx2/nX+Sy2/nY))
[1] -2.155754

et pour les degrés de liberté

> (nu=(Sx2/nX+Sy2/nY)^2/(Sx2^2/nX^2/(nX-1)+
+ Sy2^2/nY^2/(nY-1)))
[1] 36.35279

La valeur critique est obtenue en lisant dans les tables,

(car ici on a des probabilité pour un test bilatéral dans la table) comme on apprenait dans les cours de statistique au siècle passé. D’un point de vue informatique, on cherche à savoir si on est à gauche, ou à droite de la valeur critique

> qt(.05,df=nu)
[1] -1.687865

image manquante

On peut aussi calculer la p-value,

> pt(T,df=nu)
[1] 0.01889768

Si on regarde, sous R, il existe des fonctions de tests, pour comparer des moyennes. Et dans ce cas, la sortie est

> t.test(X,Y,alternative = "less")

Welch Two Sample t-test

data:  X and Y
t = -2.1558, df = 36.353, p-value = 0.0189
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
-Inf -1.772507
sample estimates:
mean of x mean of y
48.75000  56.91667

Autrement dit, on a automatiquement la p-value, et qui permet rapidement d’interpréter le test. Moralité, si on sait interpréter une p-value (et que l’on vérifié au préalable les conditions d’application d’un test), on peut faire tous les tests que l’on veut !
Si on veut faire un peu plus compliqué, on peut regarder la distribution des notes, et se demander si une loi \mathcal{N}(60,15^2) serait possible (par exemple, ça sera notre hypothèse H_0, l’hypothèse alternative étant que ce n’est pas cette loi). Pour faire ce test, il existe le test de Kolmogorov-Smirnov. La statistique de test est ici

T=\sup\{\vert \widehat{F}_n(x)-F_0(x)\vert ,x\in\mathbb{R}\}

où F_0(\cdot) est la fonction de répartition de la loi \mathcal{N}(60,15^2), et \widehat{F}_n(\cdot) est la fonction de répartition empirique

\widehat{F}_n(x)=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(x_i\leq x)

La loi de T n’est pas simple, ou moins simple qu’une loi de Student (cf Marsaglia, Tsang & Wang (2003) par exemple). En revanche, on a les p-values automatiquement,

> ks.test(Y, "pnorm", 60, 15)

One-sample Kolmogorov-Smirnov test

data:  Y
D = 0.1421, p-value = 0.5796
alternative hypothesis: two-sided

Aussi, on peut accepter ici l’hypothèse nulle. On peut d’ailleurs faire un petit dessin pour s’en convaincre,

> Femp=function(x) mean(Y<=x)
> plot(0:100,Vectorize(Femp)(0:100),type="s")
> lines(0:100,pnorm(0:100,60,15),col="red")

image manquante


Et ça va nous servir dans ce cours ? A priori oui… parce qu’on parlera du test de Student (pour tester si une variable dans une régression est significative), du test de Fisher (pour tester si plusieurs variables dans une régression sont significatives, ou plus généralement si une contrainte – linéaire – sur les coefficients peut être acceptée), du test de Chow (pour tester des ruptures dans un modèle linéaire, mais c’est un test de Fisher un peu déguisé), du test d’Anderson-Darling (pour tester si des résidus sont Gaussiens), du test de Breuch-Pagan voire le test de White (pour tester si les résidus peuvent être considérés de variance constante), du test de Durbin-Watson (pour tester s’il n’y a pas d’auto-corrélation dans la série des résidus), du test de Dickey-Fuller (pour tester si une série temporelle est – ou n’est pas – stationnaire), des tests de Franses (pour tester si une série peut être considérée comme saisonnière, ou pas), du test de Ljung-Box (pour tester si un bruit est un bruit blanc)… Et j’en oublie un paquet. Donc quand il est dit (dans le plan de cours) que le cours de statistique est un prérequis, il ne s’agit pas de l’avoir suivi, mais bel et bien de l’avoir compris, car on passera notre temps à utiliser des notions entrevues dans ce cours.