I will give, with Qichun Xu, a joint webinar for the Reinsurance Council and the Futurism Council of the Society of Actuaries, on Perspectives of Predictive Modeling with Case Studies in a few days. The slides of my talk are now available (I do recommand to open the pdf version of the slides with Acrobat, since there are animated pictures in the slides that could not be visualized below for instance). The Society of Actuaries asked specifically for a powerpoint document, so I will use screenshots of the slides for the webinar. I do encourage to open and read the pdf file for a better quality… Sorry for the inconvenience. I will upload soon lines of codes to reproduce most of the graphs. All comments and remarks are welcome.
Monthly Archives: October 2013
Halloween and candies (a ballot problem)
This year, for Halloween, a post on candies (I promise, next year I will write another post on zombies). But I don’t want to focus on the kids problems (last year, we tried to minimize their walking distance to collect as much candies as possible, with part 1 and part 2), I want to discuss my own problems. Because usually, the kids wear their costumes, and they go in the streets, they knock on the doors, while I stay at home. So I’m the one, with a bag full of candies, waiting for kids to knock on our door, and then I give them some candies (if they wear a costume). Consider the following problem. Assume that we start with red candies, and black ones, with . The thing is that no one like those black candies. What could be the probability that for the kids that will get candies after knocking at my door (with for convenience, but we will also consider the more general case where I have to many candies, , later on), the probability to get a red candy is always larger than the probability to have a black candy ? This is somehow related to the popular ballot problem, proposed (and solved) by Whitworth in 1878, but he wrote it only in the fourth edition of Choice and Chance, in 1886 (this is what the legend told us). In 1887, Joseph Bertrand proposed a similar problem, and Désiré André introduced the reflection principle to solve it. The problem is simple : consider an election between two candidates, A (who receives votes) and B (who receives votes). A wins the election (). If the ballots are cast one at a time, what is the probability that A will lead throughout the voting? For those who don’t remember the conclusion, the probability is here quite simple,
Observe that some geometry proofs were given, later on, by Aebly or Mirimanoff, both in 1923, as well as Howard Grossman in the 1950’s (see the discussion on http://academiclogbook.blogspot.ca/…). Actually, http://futilitycloset.com// produced the following geometric proof (with no clear reference),
We start at O, where no votes have been cast. Each vote for A moves us one point east and each vote for B moves us one point north until we arrive at E, the final count, (m, n). If A is to lead throughout the contest, then our path must steer consistently east of the diagonal line OD, which represents a tie score. Any path that starts by going north, through (0,1), must cut OD on its way to E.
If any path does touch OD, let it be at C. The group of such paths can be paired off as p and q, reflections of each other in the line OD that meet at C and continue on a common track to E.
This means that the total number of paths that touch OD is twice the number of paths p that start their journey to E by going north. Now, the first segment of any path might be up to m units east or up to n units north, so the proportion of paths that start by going north is n/(m + n), and twice this number is 2n/(m + n). The complementary probability — the probability of a path not touching OD — is (m –n)/(m + n).
But let’s try to solve our problem. Let and denote the number of black and red candies, respectively after the th kid git his (or her) candy. Yes, one at a time. Here, and . What we want is
Using this formulation, we recognize the ballot problem. Almost. Actually, in the original ballot problem (see Bertrand (1887)), we have to compute the probability that one candidate remains strictly ahead the other one throughout the count. With a strict condition, we get the well-known probability (given previously)
Here, ties are allowed, and we can prove (easily) that
(again, there is some nice geometric interpenetration of that result). It is also possible to get numerically that value using the following function, which will generate a trajectory, and return some indicators (with or without ties)
> red_black=function(sd){ + set.seed(sd) + vectcandy=sample(c(rep("R",r),rep("B",b))) + v1=rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")) + v2=cumsum(rev(vectcandy)=="R")<= cumsum(rev(vectcandy)=="B") + return(list(evol=cbind(rev(cumsum(rev(vectcandy)=="R")), + rev(cumsum(rev(vectcandy)=="B")),v1),list=vectcandy,test=(sum(v1)==0), + ballot=(sum(v2)==0),when=min(which(v1==1))))}
(here I compute the ballot-type index, where ties are not allowed, and the candy-type index). If we generate 100,000 scenarios, starting with 50 red and 25 black candies, we get
> r=50 > b=25 > M=sapply(1:100000,red_black) > mean(unlist(M[3,])) [1] 0.50967
which can be compared with the theoretical value
> (r+1-b)/(r+1) [1] 0.5098039
We can also get the distribution of the first time we have more black candies than red ones left (given that this event occur)
> Z=unlist(M[5,]) > Z=Z[Z<Inf] > hist(Z,breaks=seq(0,80),probability=TRUE,col="light blue", + border=NA,xlab="",main="")
There might be some analytically formula that can be derived, but I have to confess that I am becoming extremely lazy,
Assume now that this year, kids do not show up at my door (for some reason). Assume that kids show up. We can see how the probability
will change, with ,
> r=50 > b=25 > impact_n = function(n){ + red_black=function(sd,nb=n){ + set.seed(sd) + vectcandy=sample(c(rep("R",r),rep("B",b))) + v=(rev(cumsum(rev(vectcandy)=="R"))<rev(cumsum(rev(vectcandy)=="B")))[1:nb] + return(list(list=vectcandy,test=(sum(v)==0),when=min(which(v==1))))} + M=sapply(1:10000,red_black) + return(mean(unlist(M[2,])))}
Yes, not only I am too lazy to derive analytic formulas, I am so lazy that I do not try to optimize my code. Here, the evolution of the probability, as a function of is
> V=Vectorize(impact_n)(25:75) > plot(25:75,V)
Fun isn’t it? But now, I have to conclude my post, to work a little bit on my make-up : I have learnt so many thinks at the Montreal Zombie Walk a few days ago that kids willing to knock at my door will be scared to death. I guess I will keep all the candies for me this year !
More significant? so what…
Following my non-life insurance class, this morning, I had an interesting question from a student, that I will try to illustrate, and reformulate as accurately as possible. Consider a simple regression model, with one variable of interest, and one possible explanatory variable. Assume that we have two possible models, with the following output (yes, I do hide interesting parts here, but it is to get quickly to my student’s point)
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.92883 0.06391 14.534 <2e-16 *** X -0.12499 0.06108 -2.046 0.0421 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
for the first model – a GLM with some distribution, and some link function – and
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.92901 0.06270 14.817 <2e-16 *** X -0.09883 0.05816 -1.699 0.0909 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
for the second one – with another GLM, with another distribution, but the same link function (I guess I could have changed it, but it does not really matter here). Then, I got the following statement “I would like to choose the first model because the explanatory variable is more significant, and therefore, this model should have a stronger predictive power“.
That’s a nice idea, isn’t it ? Actually, I guess this is why I love teaching, because I will never be able to think about such an idea by myself. Because when you look at that statement, somehow it could make sense. Except that from my point of view, it is not valid at all. My first thought was to recall is standard example in statistical inference : you cannot not claim that a distribution is better than another one just by looking at the parameter estimates.
> fitdistr(Y,"normal") mean sd 0.93685011 0.90700830 (0.06413517) (0.04535042) > fitdistr(Y,"exponential") rate 1.06740661 (0.07547704)
Can I claim that the Gaussian distribution is better than the exponential one because parameter estimates have smaller standard deviation ? Because somehow, this is what we did when we claimed previously that the first model was better than the second one.
Let me get back on the outputs of the two regressions, and let me explain what I did. Actually, I wanted to have a story close to the one on the Gaussian versus exponential fit. So I did generate some exponential random variable,
> set.seed(5) > n=200 > U=runif(n); > Y=-log(U)
Here, we can visualize the histogram of this sample, as well as the the estimated exponential distribution
> hist(Y,proba=TRUE,col="light green",border="white",lwd=2,breaks=seq(0,5.3333333333333,by=.333333333)) > x=seq(0,6,by=.02) > lines(x,dexp(x,1/mean(Y)),col="red",lty=2)
On top of that, let us fit a gamma distribution. Using a GLM (where the regression is here on a constant – only), just to practice because later on, we will use a gamma regression on that variable
> reg0=glm(Y~1,family=Gamma(link="identity")) > a=reg0$coefficient > b=summary(reg0)$dispersion > lines(x,dgamma(x,shape=1/b,scale=a*b),col="blue")
Now, we need a covariate, to run some regressions. What I wanted is some variable slightly correlated with our previous variable. Slightly, just to make sure that our -value in the regression will be close to 5% or 10%. So here, I did generate a variable so that the pair has Clayton copula, with coefficient 0.1 (which is small, extremely small)
> a=.1 > set.seed(5) > n=200 > U=runif(n); > V=(U^(-a)*(runif(n)^(-a/(1+a))-1)+1)^(-1/a) > Y=-log(U) > X=qnorm(V)
To visualize the copula of the variables, we can use
> cop=function(u,v){ + (a+1)*(u*v)^(-(a+1))* + (u^(-a)+v^(-a)-1)^(-(2*a+1)/a) } > x=y=seq(.05,.95,by=.05) > z=outer(x,y,cop) > mat=persp(x,y,z,col="green",shade=TRUE,xlim=c(0,1),ylim=c(0,1),zlim=c(0,2),theta=-30, + ticktype ="detailed",zlab="")
We should be not far away from the independence (actually, there is a negative – significant – correlation (Pearson’s correlation)). Now, consider two models,
- a Gaussian model (here a standard linear model)
- a gamma model, with a linear link function
The outputs are the following (you will recognize the outputs given previously)
> reg1=lm(Y~X) > reg2=glm(Y~X,family=Gamma(link="identity")) > summary(reg1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.92883 0.06391 14.534 <2e-16 *** X -0.12499 0.06108 -2.046 0.0421 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.9021 on 198 degrees of freedom Multiple R-squared: 0.02071, Adjusted R-squared: 0.01576 F-statistic: 4.187 on 1 and 198 DF, p-value: 0.04206 > summary(reg2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.92901 0.06270 14.817 <2e-16 *** X -0.09883 0.05816 -1.699 0.0909 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for Gamma family taken to be 0.9086447) Null deviance: 229.72 on 199 degrees of freedom Residual deviance: 226.58 on 198 degrees of freedom AIC: 379.22 Number of Fisher Scoring iterations: 10
And here are the two predictions,
So, which model should we use? As usual, my answer will be “let’s have a look at the data” instead of looking only at tables of figures. Using some code posted a few days ago, let us visualize the two regressions. The Gaussian model is here
(for the lower part, I do not go below 0 since we do have, here, a positive variable that we would like to model) while the gamma on is here
And if we believe that the explanatory variable has no predictive power (since we can claim that the parameter is not significant in the regression), and we remove it from the regression, we get
Here, I do believe that the gamma (not to say the exponential) model is better because it is clearly more coherent with properties of the variable of interest. I trust more the confidence interval obtained above on the gamma model, than the one obtained with a Gaussian distribution. Even if the parameter in the regression is “more significant”.
Moments et variables aléatoires
Vendredi, suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur les thèmes 10, 11, 12 et 16 (tels que classifiés dans le livre de Jacques Labelle, qui sert de référence pour ce cours)
- Moments et fonction génératrice des moments #10 ACT2121-A2013-10.pdf
- Covariance et corrélation #11 ACT2121-A2013-11.pdf
- Polices d’assurance et espérances de remboursements #12 ACT2121-A2013-12.pdf
- Formule des doubles espérances #16 ACT2121-A2013-16.pdf
Réassurance
Mercredi, on finira la modélisation des coûts individuels de sinistres en évoquant la mutualisation. Si on a le temps, on parlera aussi de réassurance. Les transparents sont en ligne.
- Hurricane Katrina (US, Bahamas, Cuba, Aug. 2005), $ 72.3 billion
- Tōhoku earthquake and tsunami (Japan, March 2011), $ 35 billion
- Hurricane Andrew (US, Bahamas, August 1992), $ 25 billion
- September 11 attacks (US) $ 23.1 billion
- Northridge earthquake (US) $ 20.6 billion
- Hurricane Ike (US, Haiti, Dominican Republic, Sept. 2005) $ 20.5 billion
- Hurricane Ivan (US, Barbados, Sept. 2004) $ 14.9 billion
- Hurrican Wilman (US, Mexico, Jamaica, Oct. 2005), $ 14 billion
- Hurricane Rita (US, Cuba, Sept. 2005) $ 11.3 billion
- Hurricane Charley (US, Cuba, Jamaica) $ 9.3. billion
A titre de comparaison, les chiffres d’affaires des plus gros réassureurs (prime émise en 2010) étaient, selonhttp://www.insurancenetworking.com/…
- Munich Reinsurance Company $ 31.3 billion
- Swiss Reinsurance Company Limited $ 24.7 billion
- Hannover Rueckversicherung AG $ 15.1 billion
- Berkshire Hathaway Inc. $ 14.4 billion
- Lloyd’s $ 13 billion
- SCOR S.E. $ 8.8 billion
- Reinsurance Group of America Inc. $ 7.2 billion
- Allianz S.E. $ 5.7 billion
- PartnerRe Ltd. $ 4.9 billion
- Everest Re Group Ltd. $ 4.2 billion
De la significativité (statistique), suite
Suite à mon rapide billet sur l’article donnant les conclusions d’une étude des effets de la e-cigarette (et la comparaison avec le patch pour arrêter de fumer), j’avais entendu cette phrase “ça veut dire quoi ‘statistiquement non significatif’ ?“. Comme c’est une très bonne question, qui plus est relativement complexe, j’ai pris un peu de temps pour y répondre… Il ne s’agit pas de complexité mathématique, en s’entend car tous les cours de statistique abordent cette notion, mais plus d’une notion complexe à vulgariser. On va revenir un peu sur l’histoire des tests de significativité et sur les travaux de Ronald Fisher (pour comprendre l’origine des 5% utilisé comme seuil lors de l’utilisation de la p-value), mais il est possible d’aller voir directement les exemples.
- Un peu d’histoire
Si on veut commencer par la préhistoire, les travaux les plus anciens qui posent la question de la significativité sont ceux de John Arbuthnot, en 1710. En prenant 82 années de données de natalité, à Londres, il avait observé que chacune des années, il y avait eu plus de naissances de garçons que de filles (évoqué dans Hacking (1975)).
Il s’était demandé, avec la terminologie d’il y a plus de 300 ans, si cette différence pouvait être jugée comme ‘statistiquement significative‘. Je pense que l’on peut utiliser le terme ‘statistiquement‘ car John Arbuthnot pose la question en termes probabilistes, très clairement,
Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un ‘test de significativité‘, là encore avec notre terminologie actuelle. Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste (que l’on verra formalisé bien plus tard), en se demandant s’il est ‘extrêmement probable‘ qu’il y ait une différence entre les deux mesures. Il avait alors fait un test de comparaison des valeurs moyennes, et noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important. A la fin du XIXème siècle, Francis Edgeworth reprendra cette idée (sur la comparaison de moyennes) lorsqu’il écrira en 1885,
où on retrouve les premières idées qui serviront de base aux travaux ultérieurs. Car pour bien comprendre l’origine de la notion de significativité en statistique, il faut relire les travaux de Ronald Fisher des années 20, mais aussi de Student (William Gosset) et Karl Pearson. Ronald Fisher s’interroge, dans son papier de 1921, sur la valeur d’un coefficient de corrélation (c’est dans ce papier qu’il introduit formelle sa transformation avec une tangente hyperbolique z = {1 \over 2}\ln\left({1+r \over 1-r}\right) = \text{arctanh}(r) pour faire un test de nullité d’un coefficient de corrélation), où il note
Aussi, dans le papier, une valeur (ici une corrélation) peut être jugée comme ‘significative‘ (c’est à dire ‘significativement non nulle‘). Encore une fois, l’idée n’est pas nouvelle, Edgeworth aurait dit ‘non accidentelle‘ mais Fisher pose ici la terminologie, en parlant de ‘significativité‘. L’année suivante, en 1922, il parle de l’idée de tester la significativité, dans le cadre d’un modele de régression (on est toujours dans le contexte de la corrélation) ou en 1924, on retrouve l’idée de tester la significativité, Mais c’est surtout dans son livre Statistical Methods for Research Workers que l’on va trouver une formalisation claire (et très illustrée) de cette idée, et surtout la connecter à la notion de p-value (on ne peut pas parler de significatité sans parler de p-value). Pour l’anecdote, dans l’ouvrage de Ronald Fisher, on y parle de p-value (oui, avec une majuscule). Cela dit, la paternité de la p-value est à attribuer davantage à Pearson, par exemple dans son papier datant de 1900, ou encore, sur un autre exemple avec en plus l’interprétation usuelle en terme de significativité (sans toutefois utiliser ce terme), puisque Pearson dit qu’on a 1 chance sur 8 d’avoir une valeur aussi improbable. Plus loin, il dit aussi qu’arriver avec 1 chance sur 70 peut être qualifier de très improbable (un long débat suivra sur le seuil à partir duquel, effectivement, un événement peut être qualifié d’improbable) Mais revenons au livre de Ronald Fisher, et plus spécifiquement au chapitre 4, qui propose précisément d’étudier la ‘significativité de la différence entre deux moyennes‘
On y retrouve le test de Pearson basé sur la statistique du chi-deux, comme test d’ajustement (goodness of fit) et d’indépendance. et la comparaison de la p-value avec (ce qui va correspondre aux mythiques) 5%. Il parle plus loin d’une ‘ligne conventionnelle‘ que l’on tracera à 5%, Il faut d’ailleurs un lien entre une p-value à 5% (correspondant à un événement qui survient avec 1 chance sur 20) et le fait de s’éloigner de 2 écart-types de la moyenne pour une loi normale (qui survient avec 1 chances sur 22) Tout au long du livre, on retrouve évoqué le chiffre mythique 5% qui restera à la postérité. Le plus intéressant pour apprendre comme juger la significativité d’une valeur, on peut regarder les exemples, comme l’exemple 11 l’exemple 12, l’exemple 27, ou encore l’exemple 32 (pour en prendre juste quelques uns) (on retrouve l’utilisation de 2 écart-types pour une loi normale). La discussion amorcée dans ce livre sera poursuivie dans son papier de 1926, qui commence avec Le papier est passionnant, car il pose les bases de ce que font tous les statisticiens depuis 80 ans (il n’est pas surprenant de voir cet article reproduit dans le second volume des Breakthroughs in Statistics), où le seuil de 5% est encore utilisé. Il ne cesse tout au long de l’article de justifier cette règle de ‘un sur vingt‘, tout en autorisant d’autres valeurs, comme 2%, ou 10% (Fisher semble moins orthodoxe que certains utilisateurs de la notion de significativité) Il faudra ensuite attendre les travaux de Jerzy Neyman et Egon Pearson. qui vont axiomatiser et formaliser les tests statistiques, mais ça dépasse (de loin) l’objectif de mon billet d’aujourd’hui. Donc on va repousser à plus tard la présentation de leurs apports. Pour illustrer tout ce qu’on vient de voir, considérons deux exemples simples.
- Petit exemple, tester l’indépendance
Comme Ronald Fisher, commençons par ce qui me semble le plus simple car on n’a pas besoin de parler modèle probabiliste ici. On se demande si deux variables sont indépendantes, ou pas. Deux cas peuvent se présenter:
- on a des variables continues, et on veut voir si elles sont non-corrélées (oui, pour commencer on va simplifier, et dire que si elles sont non-corrélées… ça va se traduire par une notion que l’on peut assimiler – rapidement – à une notion d’indépendance. Promis, d’ici 6 mois il y aura des dizaines de billets sur le sujet sur le blog puisque je redonne mon cours sur la dépendance cet hiver). Je laisse les plus motivés jouer avec les données qui traînent sur le blog pour voir si la taille et le poids sont corrélées, ou utiliser un exemple que j’avais l’habitude de présenter dans mes cours de modèles de régression, où on se demande si – chez les femmes – le tour de poitrine et le tour de taille sont corrélés.
- on a des variables catégorielles, par exemple la couleur des yeux, et la couleur des cheveux. Et on se demande si les deux sont couleurs sont indépendantes, ou pas.
Regardons ce dernier exemple, avec un petit jeu de 25 observations.
> X=c(rep("Brun",8+5),rep("Blond",9+3)) > Y=c(rep("Noisette",8),rep("Bleu",5),rep("Noisette",3),rep("Bleu",9))
On peut résumer l’information dans un tableau de contingence,
> T=as.matrix(table(X,Y)) > T Y X Bleu Noisette Blond 9 3 Brun 5 8
On va essayer de comparer ce tableau avec ce qu’on devait avoir si les deux variables étaient indépendantes (car c’est ce que l’on souhaite tester : on veut savoir si le tableau que l’on a est significativement différent de ce qu’on aurait avec des variables indépendantes). Il faut ici ressortir la définition de l’indépendance entre deux variables aléatoires, \mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B) = \mathbb{P}(X\in\mathcal A) \times\mathbb{P}(Y\in\mathcal B) pour tout \mathcal A et \mathcal B. Par exemple \mathcal A peut être avoir des cheveux blonds, et \mathcal B avoir des yeux bleus. Pour rappel, on avait en tout 25 personnes ayant ces deux caractéristiques en même temps, donc empiriquement, on avait observé une probabilité de 9/25 soit 36% (c’est le terme de gauche). Pour le terme de droite, on notera que \mathbb{P}(X\in\mathcal A), c’est la probabilité d’avoir les cheveux blonds. Empiriquement, on en avait 9+3 sur 25, soit 48%. Et \mathbb{P}(Y\in\mathcal B) c’est la probabilité d’avoir les yeux bleus. Or on avait 9+5 personnes avec ces caractéristiques, soit 56%. Si les deux variables étaient indépendantes, avec la formule précédente, la probabilité d’avoir à la fois les cheveux blonds et les yeux bleus serait 48% x 56%, soit 26.88%. A comparer aux 36% observés. En fait, au lieu de comparer des pourcentages, on peut aussi comparer des nombres: comme il y a 25 personnes en tout, 26.88% signifie que dans notre groupe, 6.72 personnes devraient avoir les cheveux blonds et les yeux bleus. On peut s’amuser à faire ces calculs pour tous les \mathcal A et \mathcal B.
> Tind=chisq.test(T)$expected > Tind Y X Bleu Noisette Blond 6.72 5.28 Brun 7.28 5.72
Visuellement, on peut mettre ça en forme,
(sur le graphique de droite, les traits noirs correspondent aux nombres sous hypothèse d’indépendance, et en couleur, les valeurs observées). Bon, maintenant, il faut réfléchir un peu… On se demande si nos deux tableaux sont proches ou pas. Ou, pour utiliser notre mot de la journée, significativement proches. Alors là, ça va devenir technique le temps d’un rapide paragraphe. Dans chaque case du tableau, on compte des gens. Sur nos n personnes, on peut s’attendre à en observer un nombre aléatoire, noté N_{\mathcal{A},\mathcal{B}}, qui suit une loi binomiale N_{\mathcal{A},\mathcal{B}}\sim\mathcal{B}(n,\mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B)) (c’est la définition de la loi binomiale). On notera que quand n est grand, la loi normale est alors proche de la loi Gaussienne. Dans chaque case, on a des lois normales (non indépendantes, car les nombres sont contraints, par exemple au total, il faut n personnes), et en prenant les carrés des lois normales (centrées et réduites), on va obtenir une loi du chi-deux (je renvoie vers un précédant billet pour les aspects techniques). Prendre le carré est intéressant, et naturel, car ça fait penser à une distance Euclidienne (usuelle). On va alors poser
Q=\sum_{\mathcal{A},\mathcal{B}} \frac{[n_{\mathcal{A},\mathcal{B}}-n_{\mathcal{A},\mathcal{B}}^\perp]^2}{n_{\mathcal{A},\mathcal{B}}^\perp}
n_{\mathcal{A},\mathcal{B}} est le nombre de personnes observées, de caractéristiques (jointes) \mathcal A et \mathcal B. et n^\perp_{\mathcal{A},\mathcal{B}} est ce qu’on aurait si les variables étaient indépendantes (calculé auparavant), soit
avec les notations naturelles. Ici,
> sum((T-Tind)^2/Tind)
[1] 3.380994
Bon… on en fait quoi maintenant de ce 3.38 ? C’est là qu’on utilise le petit résultat évoqué auparavant, qui nous dit que, si n est grand, et que les variables sont effectivement indépendantes, alors Q va suivre une loi du chi-deux (en l’occurrence à 1 degré de liberté, pour tenir compte des diverses contraintes). On va alors se demander s’il est possible, ou vraisemblable, d’avoir 3.38 avec une loi du chi-deux.
> pchisq(3.380994,df=1) [1] 0.9340477 > 1-pchisq(3.380994,df=1) [1] 0.06595227
Aussi, une loi du chi-deux à un degré de liberté a 93.4% chances d’être inférieure à la valeur obtenue, ou encore
\mathbb{P}(Q>3.38) \approx 6.59\% avec Q\sim \chi^2(1)
C’est cette grandeur que l’on appelle p-value et que l’on va essayer de comparer au seuil désormais mythique de 5%,
- si cette probabilité est inférieure à 5%, on va rejeter l’hypothèse d’indépendance
- si cette probabilité est supérieure à 5%, on va retenir l’hypothèse d’indépendance (pour faire simple)
> chisq.test(T,correct=FALSE) Pearson's Chi-squared test data: T X-squared = 3.381, df = 1, p-value = 0.06595
Les plus attentifs devraient me faire remarquer que c’est vaseux mon histoire : tout ce que je raconte est valide si n est grand, et dire que 25 c’est grand… c’est limite ! Tout d’abord, je dirais que c’est ce qui est utilisé dans un paquet d’études en médecine (peut-être un peu plus, mais guère plus). Ensuite je dirais que, malheureusement, la plupart des résultats que l’on voit dans les cours de statistiques sont de ce genre, à savoir des résultats asymptotiques, valides seulement si n est grand (c’est pareil en économie, ou des tests asymptotiques sont parfois évoqués avec 25 années d’observations). Mais heureusement, en s’accrochant un peu, on peut s’affranchir de cette hypothèse. A condition de faire un peu de bootstrap. En fait, l’idée est très très simple. On peut utiliser des générateurs de nombres aléatoires pour générer des tableaux de contingences pour lesquels les variables sont indépendantes. Car générer des variables aléatoires indépendantes, c’est ce qui est le plus simple. En fait, ici, on va se contenter de mélanger les variables,
> s=1 > set.seed(s) > Xs=sample(X) > Ys=sample(Y)
En faisant une permutation, je m’assure d’avoir autant de personnes dont la couleur des cheveux est \mathcal A que dans la population initiale (je ne fais que permuter) et de personnes dont la couleur des yeux est \mathcal B. Par contre, comme je mélange indépendamment, je suis certain que les variables sont indépendantes,
> Ts=as.matrix(table(Xs,Ys)) > Ts Ys Xs Bleu Noisette Blond 7 5 Brun 7 6
Générer un tableau de contingence, c’est juste pour illustrer, car l’idée est de générer mille (ou plus encore) tableaux de contingence obtenus en simulant des variables indépendantes, et de voir ce que cette statistique du chi-deux donne
> Q=rep(NA,1000) > for(s in 1:1000){ + set.seed(s) + Xs=sample(X) + Ys=sample(Y) + Ts=as.matrix(table(Xs,Ys)) + Q[s]=sum((Ts-Tind)^2/Tind)}
Si on regarde les valeurs obtenues
> table(Q) Q 0.0509906759906761 0.337162837162837 1.06560106560107 1.92411754911755 298 262 200 133 3.38099400599401 4.81185481185481 6.9971694971695 9.00037462537462 60 33 9 3 11.9141275391275 14.489676989677 1 1
autrement dit, dans 29.8% des scénarios, on avait une distance (un statistique) valant 0.05, dans 26.2% des scénarios, on avait 0.33, et dans 20% des scénarios, on avait 1.06. Etc. Si on cumule on obtient
> cumsum(table(Q))/10 0.0509906759906761 0.33162837162837 1.06560106560107 1.92411754911755 3.38099400599401 4.81185481185481 6.9971694971695 9.00037462537462 11.9141275391275 14.489676989677
Pour rappel, on avait eu
> sum((T-Tind)^2/Tind) [1] 3.380994
qui peut effectivement être obtenu avec des variables indépendantes. Mais seulement dans 6% des scénarios. En fait, dans 95.3% des scénarios, on a eu une valeur inférieure ou égale à celle obtenue sur notre échantillon. Et dans 4.7% des scénarios, on a strictement dépassé cette valeur. On pourrait parler de p-value obtenue par bootstrap, ou rééchantillonage. On voit que l’on flirte ici avec le seuil des 5%… Visusellement, on a
> plot(c(0,as.numeric(names(table(Q)))),c(0, + cumsum(table(Q))/10),type="s",xlab="",ylab="") > abline(v=sum((T-Tind)^2/Tind),lty=2) > lines(seq(0,15,by=.01),pchisq(seq(0,15,by=.01),1)*100,col="red")
Ici, on est bien embêté pour conclure quoi que ce soit… Ce qui est amusant, c’est qu’avec un échantillon plus grand (trois fois plus grand par exemple), et les mêmes pourcentages, il n’y a plus aucune ambiguïté,
> T=T*3 > T Y X Bleu Noisette Blond 27 9 Brun 15 24
> Tind=chisq.test(T)$expected
> Tind Y X Bleu Noisette Blond 20.16 15.84 Brun 21.84 17.16
que l’on peut encore visualiser ci-dessous,
Ici, on obtient une p-value de l’ordre de 0.1%,
> chisq.test(T,correct=FALSE) Pearson's Chi-squared test data: T X-squared = 10.143, df = 1, p-value = 0.001449
On peut aussi faire du bootstrap si on trouve que 75 observations, c’est trop faible, mais la conclusion est la même,
- Petit exemple, sur la comparaison de moyenne de deux groupes
Cette fois, on a deux populations, et on se demande si les deux groupes sont ‘comparables‘. On peut penser à comparer le nombre de personnes guéries au sein de deux populations, une qui a eu un placebo, et l’autre un médicament. Cela dit, cet exemple simple est en fait compliqué car ‘guérie’ n’est peut-être pas une notion claire. Prenons plutôt des quantités que l’on peut mesurer. Disons la taille des gens, et demandons nous si dans un groupe, les gens sont significativement plus grands que dans l’autre. Je vais tricher un peu… je vais prendre deux populations que je sais être différentes, a priori. Un groupes d’hommes, et un groupes de femmes. On va noter \{x_1,\cdots,x_{n_X}\} et \{y_1,\cdots,y_{n_Y}\} ces deux échantillons, car il va falloir formaliser un minimum.
> Davis=read.table( + "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") > Davis[12,c(2,3)]=Davis[12,c(3,2)] > Davis=Davis[order(Davis$height),] > attach(Davis) > set.seed(1) > X=sample(height[sex=="F"],size=6) > Y=sample(height[sex=="M"],size=5)
On a ici 5 personnes dans le premier groupe et 6 personnes dans le second. Et on se demande si les tailles sont significativement différentes entre les deux groupes,
> (mx=mean(X)) [1] 165.5 > (sx=sd(X)) [1] 4.679744 > (my=mean(Y)) [1] 178 > (sy=sd(Y)) [1] 7.968689
On veut savoir si l’écart entre 165,5 cm et 178 cm est significativement différent de zéro (on pourrait aussi se demander s’il est significativement positif, mais ça serait une autre histoire). On va faire une hypothèse ici : on va supposer que nos échantillons sont obtenus comme des tirages de variables aléatoires, indépendantes, x_i=X_i(\omega) et y_j=Y_j(\omega). Le modèle qu’on va supposer est que la taille des gens suit une loi normale, conditionnellement à leur sexe,
\left\{\begin{array}{l}X_i\sim\mathcal{N}(\mu_X,\sigma_X^2) \\Y_j\sim\mathcal{N}(\mu_Y,\sigma_Y^2)\end{array}\right.Si on veut visualiser nos données, on peu utiliser
> u=seq(150,200,by=.25) > v1=dnorm(u,mx,sx) > v2=dnorm(u,my,sy) > plot(u,v1,col="red",type="l",ylim=c(0,.09),axes=FALSE) > axis(1) > lines(u,v2,col="blue") > points(X,rep(.075,length(X)),col="red") > points(Y,rep(.075,length(Y)),col="blue") > abline(v=mx,lty=2,col="red") > abline(v=my,lty=2,col="blue")
Le plus simple aurait été de supposer que les deux lois ont la même variance, mais quand on regarde les densités sur le graphique ci-dessus, il faut rester crédible. Bref, on va utiliser le test de Welch. L’idée est d’utiliser un petit exercice de probabilité. Soient \{X_1,\cdots,X_{n_X}\} et \{Y_1,\cdots,Y_{n_Y}\} des variables aléatoires, Gaussiennes (comme supposé juste avant de faire un dessin) et indépendantes. Posons alors
\overline{X}=\frac{1}{n_X}\sum_{i=1}^{n_X} X_i
et
s_X=\frac{1}{n_X-1}\sum_{i=1}^{n_X} [X_i-\overline{X}]^2
(avec les quantités similaires pour le second échantillon) alors si \mu_X=\mu_Y, et si on note
s = \sqrt{{s_X^2 \over n_X} + {s_Y^2 \over n_Y}}.
la variable aléatoire
T = {\overline{X} - \overline{Y} \over s[latex]<br />
va suivre une loi de Student (connue depuis les travaux de <a href="http://freakonometrics.hypotheses.org/1968">William Gosset</a>), avec</p>
[latex display="true"]\frac{(s_X^2/n_X + s_Y^2/n_Y)^2}{(s_X^2/n_X)^2/(n_X-1) + (s_Y^2/n_Y)^2/(n_Y-1)}
degrés de liberté. On notera qu'il ne s'agit pas ici d'un résultat asymptotique. Par contre, il repose intégralement sur une hypothèse - forte - de normalité des tailles des individus. L'idée est alors la même qu'auparavant : on va calculer cette statistique sur notre échantillon,
> nx=length(X) > ny=length(Y) > s=sqrt(sx^2/nx+sy^2/ny) > T=(mx-my)/s > T [1] -3.091371
et se demander s'il est vraisemblable d'avoir une telle valeur avec une loi de Student. Le nombre de degrés de liberté est ici
> d=(sx^2/nx+sy^2/ny)^2/((sx^2/nx)^2/(nx-1)+(sy^2/ny)^2/(ny-1)) > d [1] 6.218682
On pourrait alors calculer
\mathbb{P}(T>t) où T\sim\mathcal{S}td(d)
mais ça serait un test pour savoir si la différence est significativement positive, ce qui n'est pas ce qu'on avait dit que l'on cherchait : on veut ici significativement non-nulle.
\mathbb{P}(\vert T\vert> \vert t\vert) où T\sim\mathcal{S}td(d)
Si on veut visualiser un peu tout ça, on peut utiliser
> u=seq(-4,4,by=.01) > v=dt(u,df=d) > plot(u,v,type="l") > abline(v=T,col="red") > u=seq(-4,T,length=100) > polygon(c(u,rev(u)),c(dt(u,df=d),rep(0,100)),col="red",border=NA)
Pour calculer ces probabilités, on utilise
> pt(T,d)+(1-pt(-T,d)) [1] 0.02038577
et comme la loi de Student est symétrique par rapport à l’origine,
> pt(T,d)*2 [1] 0.02038577
Aussi, la p-value est de l'ordre de 2%, ce qui est suffisamment petit pour dire que, oui, nos moyennes sont différentes. C'est d'ailleurs ce que renvoie la commande
> t.test(X,Y) Welch Two Sample t-test data: X and Y t = -3.0914, df = 6.219, p-value = 0.02039 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -22.31039 -2.68961 sample estimates: mean of x mean of y 165.5 178.0
On avait vu tout à l'heure que l'idée de la valeur de 5% était liée au scénario qui se répète avec 1 chance sur 20 : si on pouvait générer d'autres échantillons, indépendants, alors 1 fois sur 20 on devrait atteindre une statistique de test de l'ordre de celle observée. Pour générer des échantillons vraiment indépendants, on avait utilisé du bootstrap, en rééchantillonnant dans notre population de test. Ici, on a bien plus que nos 11 observations, on en a en réalité 200. On peut se demande ce qui se serait passé si on avait constitué d'autres échantillons,
> Diff=rep(NA,1000) > for(s in 1:1000){ + set.seed(s) + X=sample(height[sex=="F"],size=6) + Y=sample(height[sex=="M"],size=5) + Diff[s]=(mean(X)-mean(Y))/(sqrt(sd(X)^2/nx+sd(Y)^2/ny))} > hist(Diff,col="light blue",border="white")
Sur aucun échantillon, la moyenne du premier groupe n'a dépassé la moyenne du second. Alors qu'on aurait pu penser en avoir une vingtaine (2% des mille scénarios générés). Ceci vient peut-être du fait que l'hypothèse de normalité pour les tailles peut être trop restrictive...
- Pour aller plus loin
Le billet aura été un peu long, mais on aura eu le temps de voir en détails d'où vient cette notion de significativité statistique, et ses liens avec la p-value. On aura même eu le temps de voir sur deux exemples l'importance des hypothèses sous-jacentes, en particulier si le test est asymptotique, ou s'il repose sur des hypothèses fortes de lois. Mais ce n'était qu'une courte mise en bouche. Pour aller plus loin, je suggère la lecture de Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle, M.V. (2008) Sound and fury: McCloskey and significance testing in economics, Liberman, M. (2013). "Significance", in 1885 and today et le passionnant Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Et sur les p-value, Goodman, S. (2002) A Dirty Dozen:Twelve P-Value Misconceptions. Je devrais aussi citer Ziliak, S. & McCloskey, D. (2008). The cult of statistical significance, mais je ne l'ai pas encore lu. Par contre, Xi'an en parlait sur son blog. Et si j'étais courageux, je tenterais une lecture bayésienne de ce billet.... promis, dès que je trouve un peu de temps !
Pricing Reinsurance Contracts
In order to illustrate the next section of the non-life insurance course, consider the following example1, inspired from http://sciencepolicy.colorado.edu/…. This is the so-called “Normalized Hurricane Damages in the United States” dataset, for the period 1900-2005, from Pielke et al. (2008). The dataset is available in xls format, so we have to spend some time to import it,
> library(gdata) > db=read.xls( + "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls", + sheet=1) trying URL 'http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls' Content type 'application/vnd.ms-excel' length 119296 bytes (116 Kb) opened URL ================================================== downloaded 116 Kb perl: warning: Setting locale failed. perl: warning: Please check that your locale settings: LANGUAGE = "fr_CA:fr", LC_ALL = (unset), LANG = "fr_CA.UTF-8" are supported and installed on your system. perl: warning: Falling back to the standard locale ("C").
The problem with excel spreadsheets is that some columns might have pre-specified format (here, losses are with a format 000,000,000 for instance)
> tail(db) Year Hurricane.Description State Category Base.Economic.Damage 202 2005 Cindy LA 1 320,000,000 203 2005 Dennis FL 3 2,230,000,000 204 2005 Katrina LA,MS 3 81,000,000,000 205 2005 Ophelia NC 1 1,600,000,000 206 2005 Rita TX 3 10,000,000,000 207 2005 Wilma FL 3 20,600,000,000 Normalized.PL05 Normalized.CL05 X X.1 202 320,000,000 320,000,000 NA NA 203 2,230,000,000 2,230,000,000 NA NA 204 81,000,000,000 81,000,000,000 NA NA 205 1,600,000,000 1,600,000,000 NA NA 206 10,000,000,000 10,000,000,000 NA NA 207 20,600,000,000 20,600,000,000 NA NA
To get data in a format we can play with, consider the following function,
> stupidcomma = function(x){ + x=as.character(x) + for(i in 1:10){x=sub(",","",as.character(x))} + return(as.numeric(x))}
and let’s convert those values into numbers,
> 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)
Here is the dataset we will use, from now on,
> tail(base) Year Hurricane.Description State Category Base.Economic.Damage 202 2005 Cindy LA 1 3.20e+08 203 2005 Dennis FL 3 2.23e+09 204 2005 Katrina LA,MS 3 8.10e+10 205 2005 Ophelia NC 1 1.60e+09 206 2005 Rita TX 3 1.00e+10 207 2005 Wilma FL 3 2.06e+10 Normalized.PL05 Normalized.CL05 202 3.20e+08 3.20e+08 203 2.23e+09 2.23e+09 204 8.10e+10 8.10e+10 205 1.60e+09 1.60e+09 206 1.00e+10 1.00e+10 207 2.06e+10 2.06e+10
We can visualize the normalized costs of hurricanes, from 1900 till 2005, with the 207 hurricanes (here the x-axis is not time, it is simply the index of the loss)
> plot(base$Normalized.PL05/1e9,type="h",ylim=c(0,155))
As usual, there are two components when computing the pure premium of an insurance contract. The number of claims (or here hurricanes) and the individual losses of each claim. We’ve seen – above – individual losses, let us focus now on the annual frequency.
> 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
On average, we experience about 2 (major) hurricanes per year,
> mean(db$counts) [1] 1.95283
In predictive modeling (here, we wish to price a reinsurance contract for, say, 2014), we need probably to take into account some possible trend in the hurricane occurrence frequency. We can consider either a linear trend,
> reg0 <- glm(counts~years,data=db,family=poisson(link="identity"), + start=lm(counts~years,data=db)$coefficients)
or an exponential one,
> reg1 <- glm(counts~years,data=db,family=poisson(link="log"))
We can plot those three predictions, and get a prediction for the number of (major) hurricanes in 2014,
> plot(years,counts,type='h',ylim=c(0,6),xlim=c(1900,2020)) > cpred1=predict(reg1,newdata=data.frame(years=1890:2030),type="response") > lines(1890:2030,cpred1,col="blue") > cpred0=predict(reg0,newdata=data.frame(years=1890:2030),type="response") > lines(1890:2030,cpred0,col="red") > abline(h=mean(db$counts),col="black") > (predictions=cbind(constant=mean(db$counts),linear= + cpred0[126],exponential=cpred1[126])) constant linear exponential 126 1.95283 3.573999 4.379822 > points(rep((1890:2030)[126],3),prediction,col=c("black","red","blue"),pch=19)
Observe that changing the model will change the pure premium: with a flat prediction, we expect less than 2 (major) hurricanes, but with the exponential trend, we expect more than 4…
This is for the expected frequency. Now, we should find a suitable model to compute the pure premium of a reinsurance treaty, with a (high) deductible, and a limited (but large) cover. As we will seen in class next week, the appropriate model is a Pareto distribution (see Hagstrœm (1925), Huyghues-Beaufond (1991) or a survey – in French – published a few years ago).
We can use Hill’s plot to estimate the tail index,
> library(evir) > hill(base$Normalized.PL05)
Clearly, costs of major hurricanes are heavy tailed.
Now, consider an insurance company, in the U.S., with 5% market share (just to illustrate). We will consider there \tilde Y_i= Y_i/20. The losses are given below. Consider a reinsurance treaty, with a deductible of 2 (billion) and a limited cover of 4 (billion),
For our Pareto model, consider only losses above 500 millions,
> threshold=.5 > (gpd.PL <- gpd(base$Normalized.PL05/1e9/20,threshold)$par.ests) xi beta 0.4424669 0.6705315
Keep in mind the 1 hurricane out of 8 reaches that level
> mean(base$Normalized.CL05/1e9/20>.5) [1] 0.1256039
Given that the loss exceeds 500 millions, we can now compute the expected value of the reinsurance contact,
To compute it we can use
> E <- function(yinf,ysup,xi,beta){ + as.numeric(integrate(function(x) (x-yinf)*dgpd(x,xi,mu=threshold,beta), + lower=yinf,upper=ysup)$value+ + (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf)) + }
[Nov 5th] there is a typo in the previous function, since the threshold should be used, here, as a parameter in the function, if you want to play with that function an see the impact of the threshold (see a more recent post on the same topic, but a different dataset)… but here, we do not change the threshold, so it is not a big deal.
Now, it is probably time to bring all the pieces together. We might expect a bit less than 2 (major) hurricanes per year,
> predictions[1] [1] 1.95283
and each hurricane has 12.5% chances to cost more than 500 million for our insurance company,
> mean(base$Normalized.PL05/1e9/20>.5) [1] 0.1256039
and given that an hurricane exceeds 500 million loss, then the expected repayment by the reinsurance company is (in millions)
> E(2,6,gpd.PL[1],gpd.PL[2])*1e3 [1] 330.9865
So the pure premium of the reinsurance contract is simply
> predictions[1]*mean(base$Normalized.PL05/1e9/20>.5)* + E(2,6,gpd.PL[1],gpd.PL[2])*1e3 [1] 81.18538
for a cover of 4 billion, in excess of 2.
1.This example will be found in the Reinsurance and Extremal Events chapter in the forthcoming Computational Actuarial Science with R, by Eric Gilleland and Mathieu Ribatet.
Proba, intra 2
Un rapide billet pour partager le sujet de l’examen intra de la semaine passée avec des éléments de correction (incluant des statistiques de réponse, comme pour le précédant examen intra). Toutes les remarques sur mes corrections sont les bienvenues
Couples de variables aléatoires
Vendredi, suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur le thème 13 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)
- Couples de variables aléatoires #13 ACT2121-A2013-13.pdf
GLM, non-linearity and heteroscedasticity
Last week in the non-life insurance course, we’ve seen the theory of the Generalized Linear Models, emphasizing the two important components
- the link function (which is actually the key component in predictive modeling)
- the distribution, or the variance function
Just to illustrate, consider my favorite dataset
lin.mod = lm(dist~speed,data=cars)
A linear model means here Y_i=\beta_0+\beta_1X_i+\varepsilon_i
where the residuals are assumed to be centered, independent, and with identical variance. If we visualize that linear regression, we usually see something like that
The idea here (in GLMs) is to assume Y\vertX=x\sim\mathcal{N}(\beta_0+\beta_1x,\sigma^2)
which will produce the same model as the one describe previously, based on some error term. That model can be visualized below,
attach(cars) n=2 X= cars$speed Y=cars$dist df=data.frame(X,Y) vX=seq(min(X)-2,max(X)+2,length=n) vY=seq(min(Y)-15,max(Y)+15,length=n) mat=persp(vX,vY,matrix(0,n,n),zlim=c(0,.1),theta=-30,ticktype ="detailed", box = FALSE) reggig=glm(Y~X,data=df,family=gaussian(link="identity")) x=seq(min(X),max(X),length=501) C=trans3d(x,predict(reggig,newdata=data.frame(X=x),type="response"),rep(0,length(x)),mat) lines(C,lwd=2) sdgig=sqrt(summary(reggig)$dispersion) x=seq(min(X),max(X),length=501) y1=qnorm(.95,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig) C=trans3d(x,y1,rep(0,length(x)),mat) lines(C,lty=2) y2=qnorm(.05,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig) C=trans3d(x,y2,rep(0,length(x)),mat) lines(C,lty=2) C=trans3d(c(x,rev(x)),c(y1,rev(y2)),rep(0,2*length(x)),mat) polygon(C,border=NA,col="yellow") C=trans3d(X,Y,rep(0,length(X)),mat) points(C,pch=19,col="red") n=8 vX=seq(min(X),max(X),length=n) mgig=predict(reggig,newdata=data.frame(X=vX)) sdgig=sqrt(summary(reggig)$dispersion) for(j in n:1){ stp=251 x=rep(vX[j],stp) y=seq(min(min(Y)-15,qnorm(.05,predict(reggig,newdata=data.frame(X=vX[j]),type="response"), sdgig)),max(Y)+15,length=stp) z0=rep(0,stp) z=dnorm(y, mgig[j], sdgig) C=trans3d(c(x,x),c(y,rev(y)),c(z,z0),mat) polygon(C,border=NA,col="light blue",density=40) C=trans3d(x,y,z0,mat) lines(C,lty=2) C=trans3d(x,y,z,mat) lines(C,col="blue")}
We do have two parts here: the linear increase of the average, \mathbb{E}(Y\vert X=x)=\beta_0+\beta_1x and the constant variance of the normal distribution \text{Var}(Y\vert X=x)=\sigma^2.
On the other hand, if we assume a Poisson regression,
poisson.reg = glm(dist~speed,data=cars,family=poisson(link="log"))
we have something like
This time, two things have changed simultaneously: our model is no longer linear, it is an exponential one \mathbb{E}(Y\vert X=x)=e^{\beta_0+\beta_1x}, and the variance is also increasing with the explanatory variable \text{Var}(Y\vert X=x)=e^{\beta_0+\beta_1x}, since with a Poisson regression,
Y\vert X=x\sim\mathcal{P}(e^{\beta_0+\beta_1x})
If we adapt the previous code, we get
The problem is that we changed two things when we introduced the Poisson regression from the linear model. So let us look at what happens when we change the two components independently. First, we can change the link function, with a Gaussian model but this time a multiplicative model (with a logarithm link function)
gaussian.reg = glm(dist~speed,data=cars,family=gaussian(link="log"))
which is still, here, an homoscedasctic model, but this time non-linear. Or we can change the link function in the Poisson regression, to get a linear model, but heteroscedastic
poisson.lin = glm(dist~speed,data=cars,family=poisson(link="identity"))
So this is basically what GLMs are about….
Modélisation des coûts individuels
Cette semaine, même si le réseau de l’UQAM est down, on va continuer le cours et finir la section sur la modélisation de la surdispersion pour la fréquence de sinistres. On devrait ensuite commencer la modélisation des coûts individuels. En particulier, on passera du temps autour de deux points,
- la distinction lognormale et gamma
- l’écrêtement des gros sinistres
Les transparents sont en ligne. Et la base des coûts est celle évoquée au second cours.
Equidistant points on a map
This morning, I had a comment on a recent post, regarding a graph I did upload on the blog, which was extracted from a paper now online (see http://hal.archives-ouvertes.fr/hal-00871883). Jo (from KUL, I guess I can share that piece of information) asked me
I was wondering whether you would want to share the R code for plotting figures 1 and 14? W.r.t. the former, the figure-in-figure is a nice touch; as to the latter, I am curious to know how you translated distance in km to the size parameters of the graph (par(“usr”)) for plotting the corresponding concentric circles (and the arrow indicating the radius) on top of your map.
At first, I thought I made a mistake in my plot. I mean, each time I have a question, I start to be suspicious, and I start to wonder if what I did was valid, or not. Here was the graph
Let’s make it clear: I do not draw circles here. So yes, I believe that what I did is valid. What I did is simple. First, I get the background map,
library(maps) map("world",xlim=c(130,150),ylim=c(25,45),fill=TRUE,col="light green")
Then, I need some function to compute distance from coordinates. The functions I use are
deg2rad = function(deg) return(deg*pi/180) DISTANCEDEG = function(long1, lat1, long2, lat2) { R=6371; d=acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R return(d) }
The center here will be Tokyo (東京),
X=139+45/60 Y=35+40/60
The idea now is simple: I generate a grid (here 501×501)
VX=seq(X-10,X+10,length=501) VY=seq(Y-10,Y+10,length=501) VtX=rep(VX,each=length(VY)) VtY=rep(VY,length(VX)) ZDeg=deg2rad(cbind(VtX,VtY))
I compute the distance from all those points to Tokyo, and check is the distance is larger or smaller than a given value,
L=500 D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L
If the distance is smaller than 500km, then I put a blue dot on the graph,
points(VtX[D1],VtY[D1],pch=19,cex=.2,col="light blue")
Then I use the same procedure for 250km (obviously, it is more convenient to start from larger and to go to smaller distances)
L=250 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtX[D],VtY[D],pch=19,cex=.2,col="light yellow")
Then, I did draw an arrow to ilustrate the largest distance
k=which.max(VtX[D1]) arrows((VtX[D1])[k],(VtY[D1])[k],X,Y,code=3,length=.1) text(((VtX[D1])[k]+X)/2,Y+.35,"500 km")
And now, I have the graph.
Now, the point is that it should depend on the kind of projection we use, right? So here is a function that can be used for different kind of projections (some slight changes are necessary, since the map is now centered on some point, and we cannot use standard coordinates)
library(mapproj) mapjapan = function(pr="conic",pm=45){ map("world","japan",fill=TRUE,col="light green",projection=pr, par=pm) MP=mapproject(data.frame(x=X,y=Y),projection="") Xp=MP$x Yp=MP$y VX=seq(X-10,X+10,length=501) VY=seq(Y-10,Y+10,length=501) VtX=rep(VX,each=length(VY)) VtY=rep(VY,length(VX)) MP=mapproject(data.frame(x=VtX,y=VtY),projection="") VtXp=MP$x VtYp=MP$y ZDeg=deg2rad(cbind(VtX,VtY)) L=500 D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue") L=250 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow") L=100 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue") L=50 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow") points(Xp,Yp,pch=19,cex=.4,col="red") map("world","japan",projection=pr, par=pm,add=TRUE) }
The default function here produces a map based on a conic projection,
mapjapan()
But we can also use a Bonne projection (a pseudo-conic one, named after Rigobert Bonne)
mapjapan("bonne")
mappjapan("lagrange",NULL)
and (as a last one), Albers projections,
mapjapan("albers",c(30,40))
Of course, much more projections are possible !
We do not see much here, right ? So let us play with a larger country to visualize something. Like Canada. And the distance to, say, Winnipeg.
The first thing to do is to define the coordinates of Winnipeg,
X=-(97+08/60) Y=(49+53/60)
Then, we slightly change our function
mapcanada = function(pr="conic",pm=45){ map("world","canada",fill=TRUE,col="light green",projection=pr, par=pm) MP=mapproject(data.frame(x=X,y=Y),projection="") Xp=MP$x Yp=MP$y VX=seq(X-30,X+30,length=501) VY=seq(Y-30,Y+30,length=501) VtX=rep(VX,each=length(VY)) VtY=rep(VY,length(VX)) MP=mapproject(data.frame(x=VtX,y=VtY),projection="") VtXp=MP$x VtYp=MP$y ZDeg=deg2rad(cbind(VtX,VtY)) L=2000 D1=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D1],VtYp[D1],pch=19,cex=.2,col="light blue") L=1000 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow") L=500 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light blue") L=200 D=DISTANCEDEG(ZDeg[,1],ZDeg[,2],deg2rad(X),deg2rad(Y))<L points(VtXp[D],VtYp[D],pch=19,cex=.2,col="light yellow") points(Xp,Yp,pch=19,cex=.4,col="red") map("world","canada",projection=pr, par=pm,add=TRUE) }
Now, we can have some fun
mapcanada() mapcanada("bonne",45) mapcanada("albers",c(30,40)) mapcanada("lagrange",NULL)
Fun, isn’t it ? Changing the projection will change the shape of equidistant curves.
Borders, Walls and Security Conference
This Friday, I will be at the Borders, Walls and Security Conference , chairing the session before lunch time. The conference is organized by the Center for Geopolitical Studies of the Raoul Dandurand Chair. The program in online.
Generating your own normal distribution table
It might sounds incredibly old fashion, but for my the exam for the ACT2121 probability course (to prepare for the exam P of the Society of Actuaries), I will provide a standard normal distribution table. The problem is that it is never the one we’re looking for (sometimes it is the survival function, sometimes it is the cumulative distribution function, sometimes we consider only positive values, etc). Here is the one that will be given for the exam, this Friday.
Now, here is the code to generate it.
I did use the following code to generate the table (in a latex format),
> u=seq(0,3.09,by=0.01) > p=pnorm(u) > m=matrix(p,ncol=10,byrow=TRUE
We have here the table that we wish to have in our table,
> options(digits=4) > m [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0.5000 0.5040 0.5080 0.5120 0.5160 0.5199 0.5239 0.5279 0.5319 0.5359 [2,] 0.5398 0.5438 0.5478 0.5517 0.5557 0.5596 0.5636 0.5675 0.5714 0.5753 [3,] 0.5793 0.5832 0.5871 0.5910 0.5948 0.5987 0.6026 0.6064 0.6103 0.6141 [4,] 0.6179 0.6217 0.6255 0.6293 0.6331 0.6368 0.6406 0.6443 0.6480 0.6517 [5,] 0.6554 0.6591 0.6628 0.6664 0.6700 0.6736 0.6772 0.6808 0.6844 0.6879 [6,] 0.6915 0.6950 0.6985 0.7019 0.7054 0.7088 0.7123 0.7157 0.7190 0.7224 [7,] 0.7257 0.7291 0.7324 0.7357 0.7389 0.7422 0.7454 0.7486 0.7517 0.7549 [8,] 0.7580 0.7611 0.7642 0.7673 0.7704 0.7734 0.7764 0.7794 0.7823 0.7852 [9,] 0.7881 0.7910 0.7939 0.7967 0.7995 0.8023 0.8051 0.8078 0.8106 0.8133 [10,] 0.8159 0.8186 0.8212 0.8238 0.8264 0.8289 0.8315 0.8340 0.8365 0.8389 [11,] 0.8413 0.8438 0.8461 0.8485 0.8508 0.8531 0.8554 0.8577 0.8599 0.8621 [12,] 0.8643 0.8665 0.8686 0.8708 0.8729 0.8749 0.8770 0.8790 0.8810 0.8830 [13,] 0.8849 0.8869 0.8888 0.8907 0.8925 0.8944 0.8962 0.8980 0.8997 0.9015 [14,] 0.9032 0.9049 0.9066 0.9082 0.9099 0.9115 0.9131 0.9147 0.9162 0.9177 [15,] 0.9192 0.9207 0.9222 0.9236 0.9251 0.9265 0.9279 0.9292 0.9306 0.9319 [16,] 0.9332 0.9345 0.9357 0.9370 0.9382 0.9394 0.9406 0.9418 0.9429 0.9441 [17,] 0.9452 0.9463 0.9474 0.9484 0.9495 0.9505 0.9515 0.9525 0.9535 0.9545 [18,] 0.9554 0.9564 0.9573 0.9582 0.9591 0.9599 0.9608 0.9616 0.9625 0.9633 [19,] 0.9641 0.9649 0.9656 0.9664 0.9671 0.9678 0.9686 0.9693 0.9699 0.9706 [20,] 0.9713 0.9719 0.9726 0.9732 0.9738 0.9744 0.9750 0.9756 0.9761 0.9767 [21,] 0.9772 0.9778 0.9783 0.9788 0.9793 0.9798 0.9803 0.9808 0.9812 0.9817 [22,] 0.9821 0.9826 0.9830 0.9834 0.9838 0.9842 0.9846 0.9850 0.9854 0.9857 [23,] 0.9861 0.9864 0.9868 0.9871 0.9875 0.9878 0.9881 0.9884 0.9887 0.9890 [24,] 0.9893 0.9896 0.9898 0.9901 0.9904 0.9906 0.9909 0.9911 0.9913 0.9916 [25,] 0.9918 0.9920 0.9922 0.9925 0.9927 0.9929 0.9931 0.9932 0.9934 0.9936 [26,] 0.9938 0.9940 0.9941 0.9943 0.9945 0.9946 0.9948 0.9949 0.9951 0.9952 [27,] 0.9953 0.9955 0.9956 0.9957 0.9959 0.9960 0.9961 0.9962 0.9963 0.9964 [28,] 0.9965 0.9966 0.9967 0.9968 0.9969 0.9970 0.9971 0.9972 0.9973 0.9974 [29,] 0.9974 0.9975 0.9976 0.9977 0.9977 0.9978 0.9979 0.9979 0.9980 0.9981 [30,] 0.9981 0.9982 0.9982 0.9983 0.9984 0.9984 0.9985 0.9985 0.9986 0.9986 [31,] 0.9987 0.9987 0.9987 0.9988 0.9988 0.9989 0.9989 0.9989 0.9990 0.9990 > rownames(m)=seq(0,3,b=.1) > colnames(m)=seq(0,.09,by=.01)
To put it in a nice latex format, we can use
> library(xtable) > newm=xtable(m,digits=4) > print.xtable(newm, type="latex", file="nor1.tex")
We now have a simple tex file containing a table.
\begin{table}[ht] \centering \begin{tabular}{rrrrrrrrrrr} \hline & 0 & 0.001 & 0.002 & 0.003 & 0.004 & 0.005 & 0.006 & 0.007 & 0.008 & 0.009 \\ \hline 0 & 0.5000 & 0.5040 & 0.5080 & 0.5120 & 0.5160 & 0.5199 & 0.5239 & 0.5279 & 0.5319 & 0.5359 \\ 0.1 & 0.5398 & 0.5438 & 0.5478 & 0.5517 & 0.5557 & 0.5596 & 0.5636 & 0.5675 & 0.5714 & 0.5753 \\ 0.2 & 0.5793 & 0.5832 & 0.5871 & 0.5910 & 0.5948 & 0.5987 & 0.6026 & 0.6064 & 0.6103 & 0.6141 \\ 0.3 & 0.6179 & 0.6217 & 0.6255 & 0.6293 & 0.6331 & 0.6368 & 0.6406 & 0.6443 & 0.6480 & 0.6517 \\ 0.4 & 0.6554 & 0.6591 & 0.6628 & 0.6664 & 0.6700 & 0.6736 & 0.6772 & 0.6808 & 0.6844 & 0.6879 \\ 0.5 & 0.6915 & 0.6950 & 0.6985 & 0.7019 & 0.7054 & 0.7088 & 0.7123 & 0.7157 & 0.7190 & 0.7224 \\ 0.6 & 0.7257 & 0.7291 & 0.7324 & 0.7357 & 0.7389 & 0.7422 & 0.7454 & 0.7486 & 0.7517 & 0.7549 \\ 0.7 & 0.7580 & 0.7611 & 0.7642 & 0.7673 & 0.7704 & 0.7734 & 0.7764 & 0.7794 & 0.7823 & 0.7852 \\ 0.8 & 0.7881 & 0.7910 & 0.7939 & 0.7967 & 0.7995 & 0.8023 & 0.8051 & 0.8078 & 0.8106 & 0.8133 \\ 0.9 & 0.8159 & 0.8186 & 0.8212 & 0.8238 & 0.8264 & 0.8289 & 0.8315 & 0.8340 & 0.8365 & 0.8389 \\ 1 & 0.8413 & 0.8438 & 0.8461 & 0.8485 & 0.8508 & 0.8531 & 0.8554 & 0.8577 & 0.8599 & 0.8621 \\ 1.1 & 0.8643 & 0.8665 & 0.8686 & 0.8708 & 0.8729 & 0.8749 & 0.8770 & 0.8790 & 0.8810 & 0.8830 \\ 1.2 & 0.8849 & 0.8869 & 0.8888 & 0.8907 & 0.8925 & 0.8944 & 0.8962 & 0.8980 & 0.8997 & 0.9015 \\ 1.3 & 0.9032 & 0.9049 & 0.9066 & 0.9082 & 0.9099 & 0.9115 & 0.9131 & 0.9147 & 0.9162 & 0.9177 \\ 1.4 & 0.9192 & 0.9207 & 0.9222 & 0.9236 & 0.9251 & 0.9265 & 0.9279 & 0.9292 & 0.9306 & 0.9319 \\ 1.5 & 0.9332 & 0.9345 & 0.9357 & 0.9370 & 0.9382 & 0.9394 & 0.9406 & 0.9418 & 0.9429 & 0.9441 \\ 1.6 & 0.9452 & 0.9463 & 0.9474 & 0.9484 & 0.9495 & 0.9505 & 0.9515 & 0.9525 & 0.9535 & 0.9545 \\ 1.7 & 0.9554 & 0.9564 & 0.9573 & 0.9582 & 0.9591 & 0.9599 & 0.9608 & 0.9616 & 0.9625 & 0.9633 \\ 1.8 & 0.9641 & 0.9649 & 0.9656 & 0.9664 & 0.9671 & 0.9678 & 0.9686 & 0.9693 & 0.9699 & 0.9706 \\ 1.9 & 0.9713 & 0.9719 & 0.9726 & 0.9732 & 0.9738 & 0.9744 & 0.9750 & 0.9756 & 0.9761 & 0.9767 \\ 2 & 0.9772 & 0.9778 & 0.9783 & 0.9788 & 0.9793 & 0.9798 & 0.9803 & 0.9808 & 0.9812 & 0.9817 \\ 2.1 & 0.9821 & 0.9826 & 0.9830 & 0.9834 & 0.9838 & 0.9842 & 0.9846 & 0.9850 & 0.9854 & 0.9857 \\ 2.2 & 0.9861 & 0.9864 & 0.9868 & 0.9871 & 0.9875 & 0.9878 & 0.9881 & 0.9884 & 0.9887 & 0.9890 \\ 2.3 & 0.9893 & 0.9896 & 0.9898 & 0.9901 & 0.9904 & 0.9906 & 0.9909 & 0.9911 & 0.9913 & 0.9916 \\ 2.4 & 0.9918 & 0.9920 & 0.9922 & 0.9925 & 0.9927 & 0.9929 & 0.9931 & 0.9932 & 0.9934 & 0.9936 \\ 2.5 & 0.9938 & 0.9940 & 0.9941 & 0.9943 & 0.9945 & 0.9946 & 0.9948 & 0.9949 & 0.9951 & 0.9952 \\ 2.6 & 0.9953 & 0.9955 & 0.9956 & 0.9957 & 0.9959 & 0.9960 & 0.9961 & 0.9962 & 0.9963 & 0.9964 \\ 2.7 & 0.9965 & 0.9966 & 0.9967 & 0.9968 & 0.9969 & 0.9970 & 0.9971 & 0.9972 & 0.9973 & 0.9974 \\ 2.8 & 0.9974 & 0.9975 & 0.9976 & 0.9977 & 0.9977 & 0.9978 & 0.9979 & 0.9979 & 0.9980 & 0.9981 \\ 2.9 & 0.9981 & 0.9982 & 0.9982 & 0.9983 & 0.9984 & 0.9984 & 0.9985 & 0.9985 & 0.9986 & 0.9986 \\ 3 & 0.9987 & 0.9987 & 0.9987 & 0.9988 & 0.9988 & 0.9989 & 0.9989 & 0.9989 & 0.9990 & 0.9990 \\ \hline \end{tabular} \end{table}
and the following code to get a graph, illustrating was was actually computed, in the table (see a previous post for more details)
> library("tikzDevice") > options(tikzMetricPackages = c("\\usepackage[utf8]{inputenc}", + "\\usepackage[T1]{fontenc}", "\\usetikzlibrary{calc}", "\\usepackage{amssymb}")) + tikz("normal-dist.tex", width = 8, height = 4, + standAlone = TRUE, + packages = c("\\usepackage{tikz}", + "\\usepackage[active,tightpage,psfixbb]{preview}", + "\\PreviewEnvironment{pgfpicture}", + "\\setlength\\PreviewBorder{0pt}", + "\\usepackage{amssymb}")) > u=seq(-3,3,by=.01) > plot(u,dnorm(u),type="l",axes=FALSE,xlab="",ylab="",col="white") > axis(1) > I=which((u<=1)) > polygon(c(u[I],rev(u[I])),c(dnorm(u)[I],rep(0,length(I))),col="red",border=NA) > lines(u,dnorm(u),lwd=2,col="blue") > text(-1.5, dnorm(-1.5)+.17, "$\\textcolor{blue}{X\\sim\\mathcal{N}(0,1)}$", cex = 1.5) > text(1.75, dnorm(1.75)+.25, + "$\\textcolor{red}{\\mathbb{P}(X\\leq x)=\\displaystyle{ + \\int_{-\\infty}^x \\varphi(t)dt}}$", cex = 1.5) > dev.off()
Now we have the graph in another tex file. It is possible to embed the code in a tex file, or to compile the tex file to get a pdf file. I did generate the pdf file.
Here is the tex file I finally get. It is now extremely simple to get your own normal distribution table. Now, I guess it could be possible to use sweave, or knitr. Once I’ll get a copy of Yihui’s book, I’ll try to use it to generate distribution table for my courses !
Surdispersion et comptage
Cette semaine, au cours d’assurance non-vie, on abordera la surdispersion, qui clôturera la partie du cours sur la modélisation de la fréquence de sinistres. Les transparents sont en ligne. Mais avant de parler de surdispersion, on finira la présentation des GLM. Je mets un lien vers le chapitre 15 du livre de John Fox Applied regression analysis and generalized linear models ainsi que le livre de James K. Lindsey Applying Generalized Linear Models. Je voudrais aussi renvoyer vers les notes de cours de Germán Rodríguez, avec des notes sur la régression de Poisson (avec un petit complément sur la notion de overdispersion).
Les instructions pour le second devoir seront envoyées par courriel.