Approximation de matrices

Dans le cours d’analyses des données, on va commencer l’analyse en composantes principales, en lien avec l’analyse spectrale de  ou de , ou directement l’analyse en valeurs singulières de .

Reprenons les données de l ‘enquête du quotidien L’Express,

> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",
header=TRUE)
>  M=as.matrix(BD[,2:10])
>  colnames(M)=names(BD)[2:10]
>  rownames(M)=BD$dep

Rappelons que l’on a une matrice de taille

> dim(M)
[1] 94  9

La décomposition en valeurs singulières donne des matrices

> SVD=svd(M)
> dim(M)
[1] 94  9
> dim(SVD$u)
[1] 94  9
> dim(SVD$v)
[1] 9 9

telles que . La preuve, créons la matrice produit

> M1 = SVD$u %*% diag(SVD$d) %*% t(SVD$v)

Pour vérifier que les deux matrices sont (très très) proches, on a le choix de la distance. Par exemple

à partir de la norme de Frobenius, .

> sum( (M0-M)^2 )
[1] 2.160588e-25

On peut aussi utiliser

qui est la plus grande valeur singulière de .

> max(eigen(t(M0-M)%*%(M0-M))$value)
[1] 2.018236e-25

Pour (simplement) approcher la matrice de départ par une matrice de rang plus faible, on peut  utiliser une approximation de la forme

où seules les  premières colonnes sont gardées,

> k=4
> M4 = SVD$u[,1:k] %*% diag(SVD$d[1:k]) %*% t(SVD$v[,1:k])

Les deux matrices semblent proches. Par exemple, si on regarde la première ligne on a

> M[1,1:7]
infra  vols   eco  crim   vma   vvi  camb
44.11 29.51  4.14  2.79  0.08  0.27  5.91 
> M4[1,1:6]
[1] 44.07698034 29.57490613  4.18441060   
     2.86142461  0.06670419  0.50994705

Si on regarde la distance entre ces deux matrices, on est effectivement relativement faible

> sum( (M4-M)^2 )
[1] 240.5894
> max(eigen(t(M4-M)%*%(M4-M))$value)
[1] 130.3316

cette dernière étant à comparer avec

> max(eigen(t(M)%*%M)$value)
[1] 447176.6

Cette matrice est effectivement de rang plus faible (de range 4)

> svd(M4)$d[1:8]
[1] 6.687127e+02 3.520604e+01 1.860897e+01 1.328171e+01 3.017803e-14 9.930232e-15 8.511710e-15 5.410929e-15

On peut aussi regarder une approximation de rang 2,

> k=2
> M2 = SVD$u[,1:k] %*% diag(SVD$d[1:k]) %*% t(SVD$v[,1:k])
> svd(M2)$d
[1] 6.687127e+02 3.520604e+01 1.012972e-13 9.873132e-15 6.152361e-15 4.691758e-15 3.289106e-15 1.206215e-15 1.248952e-16
> sum( (M2-M)^2 )
[1] 763.2868
> max(eigen(t(M2-M)%*%(M2-M))$value)
[1] 346.2936

On continuera pendant le cours.

Ah oui, et pour avoir davantage d’illustration, on pourra utiliser une autre base, afin de travailler sur des rangs (et non pas des valeurs réelles)

> BD = read.table("http://freakonometrics.free.fr/acp-villes.csv",sep=";",header=TRUE)
>  M=as.matrix(BD[,2:12])
>  colnames(M)=names(BD)[2:12]
>  rownames(M)=BD$Agglo

 

Collectionner les Figurines Lego

 

Tous ceux qui ont des enfants – ou qui jouent encore au Lego – le savent : régulièrement, Lego sort des séries de figurines, vendus en sachets individuels hermétiques, un peu comme les anciennes cartes à jouer. On est à la saison 15, et il y a 16 figurines par séries,

D’ordinaire, les enfants réclament un paquet, quand on va dans un magasin de jeux (genre pour acheter un cadeau pour un anniversaire) mais là, ma fille les a réclamé tous.

Bon, après le moment où je ris dans ma barbe, j’ai quand même fini par me demander combien ça pourrait coûter si on achetait des paquets, de manière à récupérer les 16 figurines.

La bonne stratégie pour faire le calcul me semble d’ordonner les personnage par ordre d’apparition,  et de modéliser le nombre de paquets achetés 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 figurine, et la ème (avec la convention que vaut 1, i.e. à l’achat du premier paquet, on a notre première figurine), on a

On notera que, si  désigne la loi géométrique,
et à partir de là, les calculs sont simples, puisque

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

i.e.

ou encore

Le nombre moyen de paquets à acheter est alors
> sum(16/(1:16)) 
[1] 54.09166

Pour ceux qui ont oublié leur calculatrice, on peut se souvenir 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(16)*16 
[1] 44.36142

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(16)+ 0.57721)*16 
[1] 53.59678

Mais je m’écarte du sujet… à 3€, il faut débourser, en moyenne, 162€ pour avoir les 16 personnages, soit plus de 10€ (en moyenne) par personnage !

> sum(16/(1:16))*3 
[1] 162.275

On est ici au même prix que l’hotel Lego friends (oui, c’est une demande récurrente)

Par simulation, on peut confirmer cet ordre de grandeur (avec ici un million de simulation de la variable),

> rpaquet=function(i) sum(1+rgeom(16,prob=(1:16)/16))
> S=unlist(lapply(1:1e6,rpaquet))
> mean(3*S)
[1] 162.1709

mais surtout, on peut avoir des quantiles, par exemple

> median(S)*3
[1] 150

après avoir déboursé 150 euros, il y a 50% de chances que l’on ait déjà toutes les figurines.

> quantile(S,.9)*3
90% 
237

 

Et il y a 1 chance sur 10 pour qu’on dépasse les 237€. Mine de rien, pour 16 petites figurines, on dépasse là le prix du Falcon Millenium,

Je pense que les enfants ne sont pas prêt que je leur rachète des petites figurines !

Trafic Journalier sur les Bus Parisiens

Sur http://opendata.stif.info/, on peut avoir accès à des trafics journaliers sur des lignes de bus en région parisienne.

> base1=read.csv(
"/home/charpentier/Téléchargements/validations-sur-le-reseau-de-surface-nombre-de-validations-par-jour-1er-sem (1).csv", sep=";")
> base1=base1[-which(base1$LIBELLE_LIGNE%in%
c("Inconnu","NON DEFINI")),]
> dim(base1)
[1] 991842      8

Il faut faire un peu attention avec le trafic, qui est ici un facteur (à cause de réponses du type ‘moins de 5’)

> base1$NB_VALD=as.numeric(
+ as.character(base1$NB_VALD))

On peut considerer une ligne de bus particulière, par  exemple entre la porte de St-Cloud et la Bibliothèque Nationale,

> id="PARIS-16 (Porte de St-Cloud)  - PARIS-13 (Bibliothèque F. Mitterrand)"

ce qui pourrait être la ligne 62,

On peut spécifier ensuite la catégorie du titre de transport, comme une carte navigo, ou une carte imagine’r (étudiants et scolaires)

> B_navigo=base1[(base1$LIBELLE_LIGNE==id)&(base1$CATEGORIE_TITRE=="NAVIGO"),]
> B_ir=base1[(base1$LIBELLE_LIGNE==id)&(base1$CATEGORIE_TITRE=="IMAGINE R"),]

Comme on va visualiser une série temporelle, on doit trier en fonction de la date

> B_navigo$date=as.Date(B_navigo$JOUR,"%d/%m/%Y")
> B_navigo=B_navigo[order(B_navigo$date),]

Sur les 6 premiers mois de 2015, on observe la série temporelle suivante (avec une carte navigo)

> plot(B_navigo$date,B_navigo$NB_VALD,type="l")

Si on regarde en fonction du jour de la semaine,  du dimanche au samedi,

> jour=as.POSIXlt(B_navigo$date)$wday
> x=B_navigo$NB_VALD
> plot(jour,x,axes=FALSE,ylim=c(-1000,max(x)))
> Z_m=tapply(as.numeric(x),as.factor(jour),mean)
> lines(as.numeric(names(Z_m)),Z_m,col="red")
> axis(2)
> text(0:6,rep(-1000,7),
+ c("D","L","M","M","J","V","S"))

on obtient

avec une tendance qui ne présente pas trop de surprise : beaucoup de monde les jours de semaine, moins le samedi, et encore moins le dimanche. Pour la carte étudiante, on a une tendance un peu différente

> B_ir$date=as.Date(B_ir$JOUR,"%d/%m/%Y")
> B_ir=B_ir[order(B_ir$date),]
> plot(B_ir$date,B_ir$NB_VALD,type="l")

avec sur les premiers mois de 2016,

On retrouve les deux périodes de vacances, et une étrange tendance: une hausse jusqu’en avril, puis une décroissance (sur les jours de semaine). Si on regarde en fonction du jour de la semaine

> jour=as.POSIXlt(B_ir$date)$wday
> x=B_ir$NB_VALD
> plot(jour,x,axes=FALSE,ylim=c(-1000,max(x)))
> Z_m=tapply(as.numeric(x),as.factor(jour),mean)
> lines(as.numeric(names(Z_m)),Z_m,col="red")
> axis(2)
> text(0:6,rep(-1000,7),
+ c("D","L","M","M","J","V","S"))

On a ici, là encore, beaucoup moins de trafic le dimanche que les jours de la semaine,

On verra – si on a le temps – pour utiliser ces données dans le cadre du cours de séries temporelles.

La vie d’un sinistre…

Pour illustrer un peu le cours d’actuariat de l’assurance dommage, un petit exemple de ce que peut être une base de sinistres (je serais ravi d’en avoir que je puisse partager avec les étudiants, dans le cours).

Regardons par exemple un sinistre,

> X

On a les dates des opérations, avec des paiements et/ou des réestimations de la charge dossier/dossier. Pour visualiser la vie de ce sinistre, on peut utiliser

> D=as.Date(X[,2],"%d/%m/%Y")
> plot(D,X[,7]+cumsum((X[,3])),type="l")
> lines(D,cumsum((X[,3])),col="red")

Le sinistre est survenu en 2007. Quand on va constituer le triangle, les montants suivants vont être affecté à la ligne correspondant à l’année de survenance 2007,

> X$year=as.POSIXlt(D)$year+1900
> aggregate(X[,3],by=list(X$year),
+ function(x) sum(x,na.rm=TRUE))
  Group.1         x
1    1997    77.597
2    1998  9568.327
3    1999  1398.074
4    2000 16266.760
5    2001 21955.175
6    2002   106.800

correspondant aux incréments observés sur le graphique ci-dessus.

Manipulation matricielle de données

Pour travailler un peu la manipulation de matrices, dans le cadre du cours d’analyse de données, quelques lignes de code.

La matrice de données est ici

> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",header=TRUE)
> M=as.matrix(BD[,2:10])
> colnames(M)=names(BD)[2:10]
> rownames(M)=BD$dep
> head(M[,1:5])
    infra  vols  eco crim  vma
D1  44.11 29.51 4.14 2.79 0.08
D2  45.97 26.29 3.97 4.54 0.07
D3  38.83 21.90 4.61 3.33 0.03
D4  49.68 30.11 5.28 4.54 0.07
D5  47.67 27.19 4.94 3.91 0.05
D6 109.21 77.71 6.87 5.25 0.24

La décomposition spectrale de

> round(t(M)%*%M)
       infra   vols   eco  crim vma  vvi  camb  
infra 308599 191433 29179 21290 642 7042 30951 
vols  191433 119961 17859 12967 404 4427 19335 
eco    29179  17859  2946  2042  59  648  2942  
crim   21290  12967  2042  1564  42  448  2129  
vma      642    404    59    42   2   19    62   
vvi     7042   4427   648   448  19  237   665  
camb   30951  19335  2942  2129  62  665  3200  
roul   55363  34745  5089  3792 113 1217  5612 
auto   24841  15747  2214  1676  53  545  2520

(j’ai viré les deux dernières colonnes pour des raisons de présentation) donne ici

> eigen(t(M)%*%M)
$values
[1] 4.471766e+05 1.239465e+03 3.462936e+02 1.764038e+02 1.303316e+02 6.623592e+01 2.624739e+01 1.761801e+01 1.564292e-01
 
$vectors

(etc). Pour rappel, les liens entre la diagonalisation et l’analyse spectrale se voit avec l’écriture

> P=eigen(t(M)%*%M)$vectors
> P%*%diag(eigen(t(M)%*%M)$values)%*%t(P)

Maintenant, deux remarques. Tout d’abord le lien entre la décomposition spectrale de cette matrice, et la décomposition en valeurs singulières

> sqrt(eigen(t(M)%*%M)$values)
[1] 668.7126681  35.2060415  18.6089656  13.2817097  11.4162872   8.1385454   5.1232210   4.1973809   0.3955113
> svd(M)
$d
[1] 668.7126681  35.2060415  18.6089656  13.2817097  11.4162872   8.1385454   5.1232210   4.1973809   0.3955113
 
$u
             [,1]          [,2]          [,3]          [,4]         [,5]          [,6]         [,7]         [,8]          [,9]
 [1,] -0.08156446 -0.0632514738 -2.596265e-02  0.0960587227 -0.064759395 -0.0430947134 -0.024719699 -0.074347281  0.0550356487
 [2,] -0.08087892  0.0541859480 -5.777492e-02 -0.0628249751 -0.099084701 -0.1311683858  0.083484434  0.011872601 -0.0195034904

(etc) et la décomspotion spectrale de l’autre produit matriciel

> sqrt(eigen(M%*%t(M))$values)
 [1] 6.687127e+02 3.520604e+01 1.860897e+01 1.328171e+01 1.141629e+01 8.138545e+00 5.123221e+00 4.197381e+00 3.955113e-01 7.546505e-06 6.656212e-06 5.285131e-06 4.270035e-06

(etc).

Maintenant, pour comprendre un peu mieux la recherche de composantes pricinpales, considérons les deux variables

> sM=M[,c(1,3)]
> plot(sM)

 

On a tout intérêt à centrer et réduire les variables (ou changer la métrique, comme vu en cours)

> sMcr=sM
> for(j in 1:2) sMcr[,j]=(sMcr[,j]-mean(sMcr[,j]))/sd(sMcr[,j])
> plot(sMcr)

Avant de faire des projections sur un axe, introduisons deux fonctions

> prod_scal=function(x,u){
+   ps=rep(NA,nrow(x))
+   for(i in 1:nrow(x)) ps[i]=sum(x[i,]*u)
+   return(ps)
+ }

> proj=function(x,u){
+   px=x
+   for(j in 1:length(u)){
+     px[,j]=prod_scal(x,u)/sqrt(sum(u^2))*u[j]  
+   }
+   return(px)
+ }

Par exemple, si on projette sur l’axe des abscisses,

> points(proj(sMcr,c(1,0)),col="blue")

On peut ensuite chercher la direction de l’axe qui nous donnera le nuage de points qui aura la plus grande inertie

> inertie=function(x) sum(x^2)
> Theta=seq(0,3.141592657,length=101)
> V=unlist(lapply(Theta,function(theta) inertie(proj(sMcr,c(cos(theta),sin(theta))))))
> plot(Theta,V,type='l')

> (angle=optim(0,function(theta) -inertie(proj(sMcr,c(cos(theta),
sin(theta)))))$par)
[1] 0.7853516

Visuellement, on obtient

> plot(sMcr)
> psMcr=proj(sMcr,c(cos(angle),sin(angle)))
> points(psMcr,col="blue")

https://freakonometrics.hypotheses.org/files/2016/01/proj1.png

On notera que l’axe qui donne l’inertie maximale est lié aux vecteurs propres de la décomposition spectrale (celui associé à la plus grande valeur propre)

> c(cos(angle),sin(angle))
[1] 0.7071397 0.7070738
> eigen(t(sMcr)%*%sMcr)
$values
[1] 159.64663  26.35337
 
$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

On continuera demain, à manipuler des matrices de données, de faire des projections, avant de commencer l’analyse en composantes principales, à proprement parler, la semaine prochaine

Manipuler des Séries Chronologiques

Ce matin en cours, on a vu un peu comment manipuler des séries temporelles, et en particulier comme jouer avec l’échelle de temps (passer de données hebdomadaire à de données mensuelles, par exemple).

Pour commencer, travaillons un peu sur les données extraites de google.fr/trends/ par exemple sur le mot clé fièvre,

report = read.csv(
"http://freakonometrics.free.fr/report-fievre.csv",
  skip=6,header=FALSE,nrows=626)
 X=report$V2
 plot(X,type="l")

Au lieu d’avoir un vecteur, on peut créer un format de série temporelle (hebdomadaire )

 X_h=ts(X,start=c(2004,3),frequency=52)
 plot(X_h)

On va essayer d’agréger les données pour les rendre mensuelles

 date1=substr(as.character(report$V1),1,10)
 date2=substr(as.character(report$V1),14,23) 
 d1=as.Date(date1,"%Y-%m-%d")
 d2=as.Date(date2,"%Y-%m-%d")
 d1[2]-d1[1]
 Vect_J=Vect_M=Vect_X=Vect_A= NULL
 for(i in 1:length(d1)){
   jours = seq(d1[i],d2[i],length=7)
   mois = as.POSIXlt(jours)$mon+1
   annees = as.POSIXlt(jours)$year+1900
   donnees = rep(X[i],7)
   Vect_J = c(Vect_J,jours)
   Vect_M = c(Vect_M,mois)
   Vect_A = c(Vect_A,annees)
   Vect_X = c(Vect_X,donnees)
 }
 base=data.frame(J=Vect_J,M=Vect_M,
                 A=Vect_A,X=Vect_X)
 head(base)
 base$AM = base$A*100+base$M
 head(base)
 Y=tapply(base$X,as.factor(base$AM),mean)
 Z=ts(as.numeric(Y[1:144]),
 start=c(2004,1),frequency=12)
 plot(Z)

On peut superposer les deux d’ailleurs (il y a peut être un léger décalage sur la série hebdomadaire)

plot(X_h,col="blue")
lines(as.numeric(time(Z)),as.numeric(Z),
col="red",lwd=2)

Continue reading Manipuler des Séries Chronologiques

Actuariat de l’Assurance Non-Vie #11

Pour le onzième (et dernier) chapitre du cours d’actuariat de l’assurance non-vie à l’ENSAE, on va finir d’étudier les modèles de provisionnement (provisions pour sinistres à payer) avec les modèles de régression et les modèles bayésien. Je n’aurais pas le temps d’aborder cette année le provisionnement sur données individuelles. Les slides sont en ligne (la version pdf téléchargeable est comme souvent plus complète que celle sur slideshare)

ENSAE, Devoir 3

Pour le troisième (et dernier) devoir du cours d’actuariat de l’assurance dommage, je vais demander de faire un calcul de provisions à l’aide de 2 méthodes vues en cours, à partir des données du fichier triangle_ensae.ods.

Je veux un rapport, de 5 pages, proposant un montant de provisions pour les 14 années passées, avec deux méthodes différentes. Dans un second temps, il faudra utiliser une méthode de simulations pour estimer le quantile à 99.5% des paiements à venir. Le rapport donnera à la fin le montant de provision complémentaire à constituer pour que la compagnie d’assurance puisse indemniser intégralement ses assurés dans 99.5% des cas.

La dead-line est le dimanche 31 janvier.

Analyse des Données et Rappels d’Algèbre Linéaire

Cette semaine, on va commencer le cours d’Analyse des Données. Pour le premier cours, on reverra les bases d’algèbre linéaire (valeurs propres, diagonalisation, décomposition en valeurs singulières, etc). Pour ceux qui souhaitent réviser un peu, je renvoie vers le support de cours de Lara Thomas, qui contient tous les résultats nécessaires. Le reste se verra au tableau, en cours.

Histoire d’illustrer les premiers concepts, on pourra utiliser la base de données securite, extrait du “palmarès des départements: où vit-on en sécurité ?” paru dans le journal L’Express (N 2589, du 15 février 2001).

Regression with Splines: Should we care about Non-Significant Components?

Following the course of this morning, I got a very interesting question from a student of mine. The question was about having non-significant components in a splineregression.  Should we consider a model with a small number of knots and all components significant, or one with a (much) larger number of knots, and a lot of knots non-significant?

My initial intuition was to prefer the second alternative, like in autoregressive models in R. When we fit an AR(6) model, it’s not really a big deal if most coefficients are not significant (but the last one). It’s won’t affect much the forecast. So here, it might be the same. With a larger number of knots, we should be able to capture small bumps that we’ll never capture with a smaller number.

Here is what a have with a small number of knots, and cubic splines

and with a larger number of knots

In order to understand what’s going on, consider a simple model, with the two splines above, in red

> set.seed(1)
> library(splines)
> x=seq(0,1,by=.01)
> v=bs(x,10)
> x2=v[,2]
> x10=v[,10]
> set.seed(1)
> y=1+3*x2+5*x10+rnorm(length(x))/4
> y_test=1+3*x2+5*x10+rnorm(length(x))/4

Note that here I have generated two sets of data, one to train a model, and one to test it.  Here, the data looks like that

> plot(x,y)

It is based on two splines,

> lines(df$x,1+3*x2+5*x10)

If we use a spline model with 10 degrees of freedom, we get

> df=data.frame(x,y)
> reg=lm(y~bs(x,10),data=df)
> summary(reg)
 
Coefficients:
            Estimate Std. Er t value Pr(>|t|)    
(Intercept)  0.91671 0.17068   5.371 6.08e-07 ***
bs(x, 10)1   0.20485 0.32696   0.627    0.533    
bs(x, 10)2   3.15593 0.22534  14.005  < 2e-16 ***
bs(x, 10)3   0.04847 0.25075   0.193    0.847    
bs(x, 10)4   0.09373 0.21597   0.434    0.665    
bs(x, 10)5   0.11624 0.22939   0.507    0.614    
bs(x, 10)6   0.24829 0.22293   1.114    0.268    
bs(x, 10)7  -0.06825 0.23498  -0.290    0.772    
bs(x, 10)8   0.19633 0.26241   0.748    0.456    
bs(x, 10)9   0.27557 0.26976   1.022    0.310    
bs(x, 10)10  4.78134 0.24116  19.826  < 2e-16 ***

which makes sense, from what we have generated. Indeed, most of the components are not significant, but the second and the tenth. We can actually test that all those components are null (at the same time)

> A=matrix(0,8,11)
> colnames(A)=names(coefficients(reg))
> A[1,2]=A[2,4]=A[3,5]=A[4,6]=A[5,7]=
+ A[6,8]=A[7,9]=A[8,10]=1
> b=rep(0,8)
> linearHypothesis(reg, A,b)
Linear hypothesis test
 
Hypothesis:
bs(x, 10)1 = 0
bs(x, 10)3 = 0
bs(x, 10)4 = 0
bs(x, 10)5 = 0
bs(x, 10)6 = 0
bs(x, 10)7 = 0
bs(x, 10)8 = 0
bs(x, 10)9 = 0
 
Model 1: restricted model
Model 2: y ~ bs(x, 10)
 
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     98 4.8766                           
2     90 4.6196  8   0.25701 0.6259  0.754

and yes, those coefficients are not significant.

> yp10=predict(reg)
> lines(df$x,yp10,col="red")

Continue reading Regression with Splines: Should we care about Non-Significant Components?