Additional thoughts about ‘Lorenz curves’ to compare models

A few month ago, I did mention a graph, of some so-called Lorenz curves to compare regression models, see e.g. Progressive’s slides (thanks Guillaume for the reference)

The idea is simple. Consider some model for the pure premium (in insurance, it is the quantity that we like to model), i.e. the conditional expected valeur

On some dataset, we have our predictions, as well as observed quantities, . The curve are obtained simply :

  • sort the observations so that

  • based on that ordering (from high risks to low risks, based on our predictions), we plot Lorenz curve

Continue reading Additional thoughts about ‘Lorenz curves’ to compare models

“L’empereur est tout nu !”, oeuvres d’art et experts

Pendant les dernières vacances on était à Venise, où on a vu de très beaux musées. Et comme j’avais l’impression de ne pas avoir parlé musées depuis longtemps, je m’étais dit que je pourrais revenir sur une petite anecdote. D’ailleurs, comme dans mon dernier billet, après le Guggenheim de NYC, on a vu le Guggenheim de Venise. Superbe. Et on a vu les deux musées qui proposent des œuvres de la collection Pinault, le Palazzo Grassi, qui hébergeait une rétrospective Martial Raysse (artiste vraiment intéressant, très pop)

et le Punta Della Dogana. Le bâtiment est réellement incroyable, et fabuleux. Des grandes pièces, hautes, avec du bois, et des briques,

Vraiment, un musée qui vaut le détour. Mais l’expo temporaire est… déroutante….

Continue reading “L’empereur est tout nu !”, oeuvres d’art et experts

De la réussite au brevet dans cinq collèges rennais

Ce matin, on a travaillé sur des données de réussite au brevet des collèges, sur Rennes (dans le cadre du cours de statistique). On dispose de données pour 5 collèges rennais, sur une douzaine d’années (entre 2001 et 2015).

  • Analyse de le temps de l’agrégation des collèges

On avait commencé par quelque chose de simple: regarder l’évolution dans le temps du taux de réussite,

> base=
+ read.table(
+ "http://freakonometrics.free.fr/brevet_rennes.csv", header=TRUE,sep=";",dec=",")
> idx=seq(4,ncol(base),by=3)
> annees=as.numeric(substr(names(base)[idx],2,5))
> College=6
> plot(annees,base[College,idx],type="b")

Continue reading De la réussite au brevet dans cinq collèges rennais

De la fréquence des ouragans aux Etats-Unis

Pour la fin du cours de statistique, on joue enfin avec des données. Après les taux de réussite au brevet des collèges, on a travaillé un peu sur la survenance des ouragans aux Etats-Unis. La base utilisée est la suivante

> base=read.table(
"http://freakonometrics.free.fr/ouragan.csv",
sep=";",header=TRUE,dec=",")
> base=base[1:207,1:7]

Le nombre annuel d’ouragan est la série suivante

>  TB <- table(base$Year)
>  years <- as.numeric(names(TB))
>  counts <- as.numeric(TB)
>  years0=(1900:2005)[which(!(1900:2005)%in%years)]
>  db <- data.frame(years=c(years,years0),
         counts=c(counts,rep(0,length(years0))))
>  X=db$counts
>  plot(db,type="h")

On peut comparer les fréquences observées avec ce que donnerait une loi de Poisson (paramètre obtenu par maximum de vraisemblance)

>  lambda=mean(X)
>  n=length(X)
>  cbind(X,dpois(0:6,lambda)*n)
  [,1]      [,2]
0   16 15.038430
1   34 29.367500
2   27 28.674870
3   13 18.665717
4    5  9.112744
5    6  3.559128
6    5  1.158396

On peut alors faire un test d’ajustement,

>  library(vcd)
>  fit=goodfit(X,"poisson")
>  plot(fit)

>  summary(fit)
 
	 Goodness-of-fit test for poisson distribution
 
                      X^2 df   P(> X^2)
Likelihood Ratio 14.17686  5 0.01452415

On voit que le test d’ajustement d’une loi Poisson n’est pas très concluant ici. Notons que ni la loi binomiale

>  fit=goodfit(X,"binomial")
>  plot(fit)
>  summary(fit)
 
	 Goodness-of-fit test for binomial distribution
 
                      X^2 df     P(> X^2)
Likelihood Ratio 46.21196  5 8.222959e-09

ni la loi binomiale négative ne semble marcher.

>  fit=goodfit(X,"nbinomial")
>  plot(fit)
>  summary(fit)
 
	 Goodness-of-fit test for nbinomial distribution
 
                      X^2 df   P(> X^2)
Likelihood Ratio 10.91507  4 0.02753523

On peut faire un test de stabilité temporelle du paramètre d’une loi de Poisson, à l’aide du test de Przyborowski & Wilenski (1940).  Un ancien package permettait de mettre en place ce test. Mais on peut heureusement utiliser une autre fonction, avec une version différente du test (on pourra aussi consulter Krishnamoorthy & Thomson (2004) pour aller plus loin). L’idée du premier test est assez simple: notons que si sur la première période (de durée ), on a un processus de Poisson de paramètre  et sur la seconde période (de durée ), alors

et

Un test de  s’écrit  Notons que la loi
conditionnelle de  sachant  suit une loi binomiale, donc le paramètre de probabilité est fonction de ce ratio. Plus précisément, conditionnellement à,  suit un loi binomiale  avec

A partir de là, il est facile de mettre en place une stratégie de test. Une autre stratégie est d’utiliser

>  seuil=1960
>  X1=db[db$years<seuil,"counts"]
>  X2=db[db$years>=seuil,"counts"]
>  N1=length(db$years<seuil)
>  N2=length(db$years>=seuil)
>  poisson.test(c(sum(X1),sum(X2)),c(N1,N2))
 
	Comparison of Poisson rates
 
data:  c(sum(X1), sum(X2)) time base: c(N1, N2)
count1 = 85, expected count1 = 103.5, p-value = 0.01216
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.5218662 0.9265666
sample estimates:
rate ratio 
 0.6967213

 

Bref, on peut accepter que le nombre annuel d’ouragan n’est pas un processus de Poisson homogène dans le temps…

Autres Données pour le cours de Statistiques

Encore quelques données pour faire un TP de statistique. Ce sont des données de coûts d’ouragans, aux Etats-Unis,

> library(gdata)
> db=read.xls(
+     "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls",
+     sheet=1)
> stupidcomma = function(x){
+  x=as.character(x)
+  for(i in 1:10){x=sub(",","",as.character(x))}
+  return(as.numeric(x))}
> base=db[,1:4]
> base$Base.Economic.Damage=
Vectorize(stupidcomma)(db$Base.Economic.Damage)
> base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05)
> base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05)

Lire la base est un peu pénible, à cause de virgules qui traînent comme séparateur de milliers. Mais peu importe, le petit code permet de nettoyer la base.

  • Fréquence d’ouragans

La première étape pourra être de travailler sur la fréquence. On pourra ajuster des lois de Poisson, ou binomiales, ou binomiales négatives, et faire des tests de stabilité du paramètre, avec le temps (on pourra  couper en 4 ou 5 blocs).

> TB <- table(base$Year)
> years <- as.numeric(names(TB))
> counts <- as.numeric(TB)
> years0=(1900:2005)[which(!(1900:2005)%in%years)]
> db <- data.frame(years=c(years,years0),
+       counts=c(counts,rep(0,length(years0))))
> db[88:93,]
   years counts
88  2003      3
89  2004      6
90  2005      6
91  1902      0
92  1905      0
93  1907      0
> plot(years,counts,type='h')

  • Coût des ouragans

Dans un second temps, on pourra travailler sur le coût des ouragans

> plot(base$Normalized.PL05/1e9,type="h")

On pourra tenter d’ajuster des lois, log normale, gamma, Pareto, etc. Faire des tests d’ajustement. Et tester la stabilité des paramètres, dans le temps (en coupant en 2 blocs, par exemple, ou davantage).

Données pour le cours de Statistique

Pour les derniers cours de Statistique, j’avais envoyé  une feuille avec quelques exercices, à regarder. J’en profile pour mettre en ligne les exercices sur l’inférence, et ceux sur les tests. Je mets aussi en ligne un tableau contenant les résultats de cinq établissements sur Rennes, au Brevet des Collèges. On pourra utiliser ces données pour faire un peu de statistique.

Overview on Multivariate Distributions

In June 2016, with Olivier L’Haridon, we will organize a (small) conference, in Rennes, on risk models in a multi-attribute framework. In order to fully enjoy the workshop (more to come on the blog), we will organize every month an internal workshop on that topic. We will start tomorrow afternoon, 13:00-14:30, and I will give a brief talk on multivariate distributions, with an emphasis on spherical / elliptical distributions, distributions on the simplex, and copulas. Slides are now online,

‘Mathiness’ et Assurance

Au printemps dernier, l’économiste Paul Romer avait initié une discussion passionnante autour de la notion de “mathiness” (“mathitude” pourrait être la traduction la plus simple pour reprendre l’expression de Lavigne (2015)), posant la question de la place du formalisme mathématique en économie. Ou plus précisément, des liens entre les théorèmes mathématiques et la réalité économique. Et la principale critique faite par Paul Romer est que les mathématiques sont trop souvent utilisées pour obscurcir, plutôt qu’éclairer un débat. Cette discussion n’est pas nouvelle, mais il est intéressant de voir qu’elle survient à un moment où les assureurs s’interrogent, eux aussi, sur les modèles actuariels.

Science économique, et économie

Le point de départ de la discussion de Romer (2015) est qu’il est surprenant qu’en plus de 20 ans, la théorie de la croissance n’ait toujours pas abouti à un consensus. Son explication repose sur le concept de mathiness. Pour reprendre l’expression de Tim Harford, “economics has a reputation for producing rigourous nonsense“. Dans les années 40, Paul Samuelson avait initié un mouvement visant à clarifier les arguments des économistes, en cherchant des fondements plus « scientifiques ». La route de l’enfer serait-elle pavée de bonnes intentions ?

Continue reading ‘Mathiness’ et Assurance

Données TP ANOVA

Demain, en cours de statistique mathématique, on verra l’Analyse de la Variance, à un et deux facteurs. Les données que nous utiliserons sont des durées d’appels de téléphone cellulaire,

> download.file(
"http://freakonometrics.free.fr/celulaire.Rdata","cellulaire")
> load("cellulaire")
> head(cellulaire)
  Jour Duree  EntSor
1  LUN    30 Sortant
2  LUN    47 Entrant
3  LUN     6 Entrant
4  LUN     6 Sortant
5  LUN    17 Sortant
6  LUN   116 Entrant
> tail(cellulaire)
   Jour Duree  EntSor
60  SAM    46 Sortant
61  DIM     4 Sortant
62  DIM    44 Sortant
63  DIM    10 Sortant
64  DIM    63 Sortant
65  DIM   193 Entrant

Profile Likelihood

Consider some simulated data

> set.seed(1)
> x=exp(rnorm(100))

Assume that those data are observed i.id. random variables with distribution, with . The natural idea is to consider the maximum likelihood estimator

For instance, consider some maximum likelihood estimator,

> library(MASS)
> (F=fitdistr(x,"gamma"))
     shape       rate   
  1.4214497   0.8619969 
 (0.1822570) (0.1320717)
> F$estimate[1]+c(-1,1)*1.96*F$sd[1]
[1] 1.064226 1.778673

Here, we have an approximated (since the maximum likelihood has an asymptotic Gaussian distribution) confidence interval for . We can use numerical optimization routine to get the maximum of the log-likelihood function

> log_lik=function(theta){
+   a=theta[1]
+   b=theta[2]
+   logL=sum(log(dgamma(x,a,b)))
+   return(-logL)
+ }

> optim(c(1,1),log_lik)
$par
[1] 1.4214116 0.8620311
 
$value
[1] 146.5909

And we have the same value.

Now, what if we care only about , and not . The we can use profile likelihood. The idea is to solve

i.e.

or, equivalently,

> prof_log_lik=function(a){
+   b=(optim(1,function(z) -sum(log(dgamma(x,a,z)))))$par
+   return(-sum(log(dgamma(x,a,b))))
+ }

> vx=seq(.5,3,length=101)
> vl=-Vectorize(prof_log_lik)(vx)
> plot(vx,vl,type="l")
> optim(1,prof_log_lik)
$par
[1] 1.421094
 
$value
[1] 146.5909

A few weeks ago, we have mentioned the likelihood ratio test, i.e.

The analogous can be obtained here, since

(the 1 comes from the fact that  is a one-dimensional coefficient). The (technical) proof can be found in Suhasini Subba Rao’s notes (see also Section 4.5.2 in Antony Davison’s Statistical Models). From that property, we can easily obtain a confidence interval for 

Hence, from our sample, we get the following 95% confidence interval,

> abline(v=optim(1,prof_log_lik)$par,lty=2)
> abline(h=-optim(1,prof_log_lik)$value)
> abline(h=-optim(1,prof_log_lik)$value-qchisq(.95,1)/2)
 
> segments(F$estimate[1]-1.96*F$sd[1],
-170,F$estimate[1]+1.96*F$sd[1],-170,lwd=3,col="blue")
> borne=-optim(1,prof_log_lik)$value-qchisq(.95,1)/2
> (b1=uniroot(function(z) Vectorize(prof_log_lik)(z)+borne,c(.5,1.5))$root)
[1] 1.095726
> (b2=uniroot(function(z) Vectorize(prof_log_lik)(z)+borne,c(1.25,2.5))$root)
[1] 1.811809

that can be visualized below,

> segments(b1,-168,b2,-168,lwd=3,col="red")

In blue the obtained obtained using the asymptotic Gaussian property of the maximum likelihood estimator, and in red, the obtained obtained using the asymptotic chi-square distribution of the log (profile) likelihood ratio.

Cite this article as: Arthur Charpentier, "Profile Likelihood," in Freakonometrics, 16/11/2015, http://freakonometrics.hypotheses.org/20573.