Freakonometrics

A la suite des remarques des élèves de M2 lors de la réunion pédagogique, je rajoute un TP dans le cours d'actuariat de M1, qui aura lieu mercredi après midi. Nous ferons un peu ensemble ce que j'ai demandé pour les DM sur d'autres bases, et d'autres contrats. Pour commencer, travaillons sur le DM2, i.e. le modèle collectif. Dans un premier temps, un premier billet pour modéliser la fréquence de sinistres.
Rappelons que nous avons présenté 3 lois de base pour modéliser la fréquence,

• la loi binomiale, i.e. 
• la loi de Poisson, i.e. 
• la loi binomiale négative, i.e. 

Afin d'opérer rapidement une première sélection, on peut calculer espérance et variance empiriques,
> N=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/base-nb.txt",header=TRUE)$x
> mean(N)
[1] 0.1216735
> var(N)
[1] 0.1831262
Avant de se lancer dans toute modélisation, regardons également la distribution du nombre de sinistres par contrat
> (freq.empirique = table(N))
N
    0     1     2     3     4     5     7    16    21
35549  3008   640   119    22     2     1     1     1
En l'occurence, la loi binomiale négative semblerait plus adaptée. Mais auparavant, regardons les autres lois.Dans le cas de la loi de Poisson, l'estimateur du maximum de vraisemblance et l'estimateur par la méthode des moments coïncide. Le paramètre est ici
> lambda=mean(N)
Si l'on compare la distribution empirique avec la loi de Poisson, on obtient
> (freq.theorique = length(N)*dpois(as.numeric(names(freq.empirique)),lambda))
[1] 3.483576e+04 4.238589e+03 2.578619e+02 1.045832e+01 3.181251e-01
[6] 7.741478e-03 2.728767e-06 3.841856e-24 4.195623e-35

On peut aussi utiliser différents tests pour vérifier la qualité de l'ajustement
> library(vcd)  
Le chargement a nécessité le package : MASS
Le chargement a nécessité le package : grid
Le chargement a nécessité le package : colorspace
> gof = goodfit(N,type= "poisson",method= "ML")
> plot(gof,main="Ajustement d'une loi de Poisson")

Le test évoqué ici est celui du chi-deux, qui formellement donne les résultats suivants (en bricolant à la main pour obtenir une distribution empirique),
> x=seq(0,25)
> y=table(N)/length(N)
> freq.empirique=rep(0,26)
> freq.empirique[1:6]=y[1:6]
> freq.empirique[8]=y[7]
> freq.empirique[16+1]=y[8]
> freq.empirique[22]=y[9]
>  (freq.theorique1 = dpois(x,lambda))
 [1] 8.854374e-01 1.077343e-01 6.554202e-03 2.658242e-04 8.085939e-06
 [6] 1.967689e-07 3.990259e-09 6.935839e-11 1.054885e-12 1.426128e-14
[11] 1.735219e-16 1.919365e-18 1.946132e-20 1.821482e-22 1.583044e-24
[16] 1.284096e-26 9.765029e-29 6.989089e-31 4.724371e-33 3.025425e-35
[21] 1.840570e-37 1.066422e-39 5.897966e-42 3.120114e-44 1.581813e-46
[26] 7.698588e-49 
> chisq.test(x=freq.empirique,p=freq.theorique1)
        Chi-squared test for given probabilities
data:  freq.empirique
X-squared = 6.058094e+29, df = 25, p-value < 2.2e-16
Bref, on a du mal à accepter l'ajustement de la loi de Poisson (mais on s'y attendait un peu).
Pour la loi binomiale négative, les méthodes classiques (méthode des moments et méthode du maximum de vraisemblance) n'ont pas de solution analytique. Il faut donc faire des calculs. Pour le maximum dde vraisemblance, on peut utiliser le code suivant,
> library(MASS)
> fitdistr(N,"negative binomial")
      size           mu    
  0.287203791   0.121677345
 (0.013202351) (0.002098181)
Mais on peut aussi tenter un test de goodness of fiit,
> gof = goodfit(N,type= "nbinomial",method= "ML")
> plot(gof,main="Ajustement d'une loi Binomiale Négative")
On peut là aussi faire un test du chi-deux,
> x=seq(0,100)
> y=table(N)/length(N)
> freq.empirique=rep(0,101)
> freq.empirique[1:6]=y[1:6]
> freq.empirique[8]=y[7]
> freq.empirique[16+1]=y[8]
> freq.empirique[22]=y[9]
> parametres= fitdistr(N,"negative binomial")$estimate
>  (freq.theorique2 = dnbinom(x,parametres[1],mu=parametres[2]))
  [1] 9.035266e-01 7.722249e-02 1.479019e-02 3.355599e-03 8.206336e-04
  [6] 2.093949e-04 5.491026e-05 1.467661e-05 3.978407e-06 1.090153e-06
[...]
 [96] 1.104459e-52 3.262314e-53 9.636853e-54 2.846935e-54 8.411084e-55
[101] 2.485180e-55
> chisq.test(x=freq.empirique,p=freq.theorique2)
        Chi-squared test for given probabilities
data:  freq.empirique
X-squared = 2237.183, df = 100, p-value < 2.2e-16

autrement dit, on accepte là aussi difficilement l'hypothèse de loi binomiale négative, même si la distance du chi-deux est beaucoup plus faible qu'avec la loi de Poisson. Bref, comme il faut retenir un modèle, on acceptera le modèle suivant
dont la  moyenne théorique vaut 0,121, et la variance théorique 0,17, ce qui n'est pas trop éloigné de
> mean(N)
[1] 0.1216735
> var(N)
[1] 0.1831262