Tag Archives: Student

De la pratique de la régression

Depuis le début de la session, j’ai imposé une petite innovation, en donnant, environ une semaine sur deux, un petit exercice (obligatoire mais non noté) avant le cours, en vue de forcer à réfléchir (et de donner des éléments de réponse). Par exemple pour le premier cours, il fallait “prévoir” une valeur manquante, et le but était de montrer que, naturellement, on choisit la valeur moyenne.

Pour demain, j’avais posé un exercice un peu plus compliqué, sachant qu’on avait vu, lors du dernier cours, comme faire une régression linéaire, et qu’on avait fini en discutant les tests simples (en lien avec la significativité) et les tests multiples. Pour l’exercice, j’avais mis en ligne une petite base de données,

download.file("http://freakonometrics.free.fr/data3.RData","data3.RData")
load("data3.RData")
str(df)
'data.frame':	147 obs. of  10 variables:
 $ Y : num  11.72 15.91 14.19 11.15 8.31 ...
 $ X1: num  1.33 3.18 0.28 2.08 0.11 1.67 1.97 1.27 4.38 0.52 ...
 $ X2: num  3.66 3.75 3.32 2.68 4.97 2.98 4.56 1.78 2.83 6.36 ...
 $ X3: num  1.41 3.01 0.34 2.19 0.25 1.69 2.01 1.25 4.41 0.43 ...
 $ X4: num  -3.53 -4.46 -3.35 -7.54 -7.02 -2.53 -6.1 -5.99 -3.92 -5.84 ...
 $ X5: num  0.57 0.01 -0.7 1.62 -0.95 -1.37 1.18 -0.72 2.63 -1.63 ...
 $ X6: num  -0.82 1.2 3.03 -0.91 -1.6 1.77 1 -1.33 1.31 -0.7 ...
 $ X7: num  1.01 0.06 2.02 3.63 2.66 2.53 1.29 3.5 1.17 1.8 ...
 $ X8: num  8.31 8.52 9.78 7.34 7.26 ...
 $ X9: num  6.04 6.53 7.52 5.61 4.52 6.06 6.2 5.99 6.93 4.38 ...

Je vais mettre ici les questions que je posais, et donner des pistes de réflexions, non pas sur les réponses attendues (je n’attends rien de cet exercice à part une réflexion), mais sur la discussion que peut amener chacune des questions,

  • Faîtes un modèle linéaire pour expliquer Y en fonction des neuf variables explicatives. Combien de variables explicatives garderiez-vous ?

Commençons par faire une régression sur toutes les variables

summary(lm(Y~., data=df))
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.046133   2.193448   3.668 0.000349 ***
X1           0.342293   0.915894   0.374 0.709186    
X2          -0.040479   0.103409  -0.391 0.696073    
X3           1.683875   0.897278   1.877 0.062693 .  
X4          -0.009254   0.062382  -0.148 0.882295    
X5          -1.085367   0.113840  -9.534  < 2e-16 ***
X6           0.983207   0.111830   8.792 5.49e-15 ***
X7          -0.015646   0.087483  -0.179 0.858327    
X8           0.012165   0.094756   0.128 0.898033    
X9           0.172210   0.180605   0.954 0.342005    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.019 on 137 degrees of freedom
Multiple R-squared:  0.9228,	Adjusted R-squared:  0.9178 
F-statistic:   182 on 9 and 137 DF,  p-value: < 2.2e-16

On a, dans la sortie, une dizaine de tests de Student qui sont évoqués, correspondant au test de l’hypothèse H_0:\beta_j=0 (contre l’hypothèse alternative (bilatérale) H_1:\beta_j\neq 0), dans un modèle de la forme y_i=\beta_0+\beta_1x_{1,i}+\dots+\beta_9x_{9,i}+\varepsilon_iC’est ce qu’on appelle le test de significativité de la variable x_j (oui, comme on l’a vu en cours, on peut dire le test parce que les autres tests classiques – Fisher, ou Wald – sont équivalent – sauf qu’au lieu de regarder t on regarde la statistique t^2 – ce qui présente l’avantage de voir la statistique de test comme une forme de distance à l’hypothèse H_0 : si c’est trop grand, on rejette…). Avec un seuil d’acceptation de l’ordre de 5%, on nous dit que 6 variables ne sont pas significatives. Mais gardons bien en mémoire que le test de significativité de x_j est fait ici en supposant que toutes les autres variables restent dans le modèle. Autre chose: avec un seuil d’acceptation de l’ordre de 10%, une des variables (la troisième) peut être vue comme significative.

Faisons un test multiple, pour savoir si on peut supprimer 6 des 9 variables explicatives possibles (faisons le à la main, inutile d’aller chercher un package pour le faire)

reg1 = lm(formula = Y ~ ., data = df)
reg0 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg0)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.11031    0.16788   54.27   <2e-16 ***
X3           1.96784    0.07604   25.88   <2e-16 ***
X5          -0.99296    0.08751  -11.35   <2e-16 ***
X6           1.08761    0.05934   18.33   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: < 2.2e-16 anova(reg0,reg1) Analysis of Variance Table Model 1: Y ~ X3 + X5 + X6 Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 Res.Df RSS Df Sum of Sq F Pr(>F)
1    143 146.16                           
2    137 142.12  6    4.0338 0.6481 0.6916

Le test de Fisher nous dit qu’on peut accepter l’hypothèse que les 6 coefficients sont nuls – ici H_0:\beta_1=\beta_2=\beta_4=\beta_7=\beta_8=\beta_9=0(qui est un test multiple, contrairement au test de Student précédant). Mais il nous dit aussi, qu’individuellement, les trois variables restants semblent significatives. Donc j’aurais tendance à garder 3 variables explicatives (je ne parle pas de la constante : on garde toujours la constante – qui n’explique pas grand chose, sauf la valeur moyenne de y).

  • faites une prévision pour un individu dont on sait que X3=1, X5=1 et X6=8

(en réalité, la vraie question que j’ai posée contenait une typo ce qui la rendait vicieuse parce que ce n’est pas le modèle qu’on vient de construire… mais pour commencer, regardons cette question)

On vient de calibrer ce modèle, donc il suffit de faire une prévision,

predict(reg0, newdata=data.frame(X3=1, X5=1, X6=8))
       1 
18.78603

Mais ce n’était pas la vraie question…

  • faites une prévision pour un individu dont on sait que X1=3, X5=1 et X6=8

Je pense que pour répondre à cette question, il convient d’oublier tout ce qu’on vient de voir. On nous donne trois informations, et on va voir si on peut les exploiter. Autrement dit, on va commencer par regarder la régression sur ces trois variables,

reg156 = lm(formula = Y ~ X1+X5+X6, data = df)
summary(reg156)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.06688    0.17222   52.65   <2e-16 ***
X1           1.98342    0.07800   25.43   <2e-16 ***
X5          -1.01438    0.08947  -11.34   <2e-16 ***
X6           1.07824    0.06039   17.86   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.026 on 143 degrees of freedom
Multiple R-squared:  0.9183,	Adjusted R-squared:  0.9166 
F-statistic: 535.8 on 3 and 143 DF,  p-value: < 2.2e-16

Le modèle est ici bon, les trois variables étant significatives. Cela dit, je dis qu’il est “bon” mais on verra vendredi comment discuter davantage ce point… En tous cas, on peut tenter de faire une prévision, et on obtient

predict(reg156, newdata=data.frame(X1=3, X5=1, X6=8))
       1 
22.62867
  • faites une prévision pour un individu dont on sait que X1=3, X3=2, X5=1 et X6=8

Comme auparavant, on nous donne 4 informations… on va regarder le modèle avec les 4 variables

reg1356 = lm(formula = Y ~ X1+X3+X5+X6, data = df)
summary(reg1356)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.10458    0.17131  53.148   <2e-16 ***
X1           0.16362    0.89028   0.184    0.854    
X3           1.80663    0.88052   2.052    0.042 *  
X5          -0.99558    0.08895 -11.192   <2e-16 ***
X6           1.08649    0.05986  18.152   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.014 on 142 degrees of freedom
Multiple R-squared:  0.9207,	Adjusted R-squared:  0.9184 
F-statistic: 411.9 on 4 and 142 DF,  p-value: < 2.2e-16

Cette fois, une des variables n’est pas significative (la première) donc on devrait l’enlever: on fait alors la régression juste sur les trois autres variables

reg356 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg356)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.11031    0.16788   54.27   <2e-16 ***
X3           1.96784    0.07604   25.88   <2e-16 ***
X5          -0.99296    0.08751  -11.35   <2e-16 ***
X6           1.08761    0.05934   18.33   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: < 2.2e-16

Cette fois, le modèle est “bon” et on peut alors faire la prévision

predict(reg356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.75388

A titre d’information, si on avait fait la prévision sans enlever la variable non significative, on aurait obtenu la valeur suivante

predict(reg1356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.90501

qui est globalement assez proche. On aura l’occasion d’en reparler en cours, quand on se demandera s’il est plus grave d’avoir une valeur non-significative dans la régression, ou d’oublier une variable importante… Mais encore une fois, le but de ces petits exercices est d’appliquer ce qu’on a vu en cours, et d’introduire des questions auxquelles j’apporterai des réponses au prochain cours !

Probabilistic Foundations of Econometrics, part 4

This post is the fourth one of our series on the history and foundations of econometric and machine learning models. Part 3 is online here.

Goodness of Fit, and Model

In the Gaussian linear model, the determination coefficient – noted R^2 – is often used as a measure of fit quality. It is based on the variance decomposition formula \underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{total variance}}=\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{\text{residual variance}}+\underbrace{\frac{1}{n}\sum_{i=1}^n (\widehat{y}_i-\bar{y})^2}_{\text{explained variance}} The R^2 is defined as the ratio of explained variance and total variance, another interpretation of the coefficient that we had introduced from the geometry of the least squares R^2= \frac{\sum_{i=1}^n (y_i-\bar{y})^2-\sum_{i=1}^n (y_i-\widehat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2}The sums of the error squares in this writing can be rewritten as a log-likelihood. However, it should be remembered that, up to one additive constant (obtained with a saturated model) in generalized linear models, deviance is defined by {Deviance}(\widehat{\beta}) = -2\log[\mathcal{L}] which can also be noted Deviance(\widehat{\mathbf{y}}). A null deviance can be defined as the one obtained without using the explanatory variables \mathbf{x}, so that \widehat{y}_i=\overline{y}. It is then possible to define, in a more general context (with a non-Gaussian distribution for y)R^2=\frac{{Deviance}(\overline{y})-{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}=1-\frac{{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}However, this measure cannot be used to choose a model, if one wishes to have a relatively simple model in the end, because it increases artificially with the addition of explanatory variables without significant effect. We will then tend to prefer the adjusted R^2,\bar R^2 = {1-(1-R^{2})\cdot{n-1 \over n-p}} = R^{2}-\underbrace{(1-R^{2})\cdot{p-1 \over n-p}}_{\text{penalty}}where p is the number of parameters of the model. Measuring the quality of fit will penalize overly complex models.

This idea will be found in the Akaike criterion, where AIC=Deviance+2\cdot p or in the Schwarz criterion, BIC=Deviance+log(n)\cdot p. In large dimensions (typically p>\sqrt{n}), we will tend to use a corrected AIC, defined by AIC_c=Deviance+2⋅p⋅n/(n-p-1) .

These criterias are used in so-called “stepwise” methods, introducing the set methods. In the “forward” method, we start by regressing to the constant, then we add one variable at a time, retaining the one that lowers the AIC criterion the most, until adding a variable increases the AIC criterion of the model. In the “backward” method, we start by regressing on all variables, then we remove one variable at a time, removing the one that lowers the AIC criterion the most, until removing a variable increases the AIC criterion from the model.

Another justification for this notion of penalty (we will come back to this idea in machine learning) can be the following. Let us consider an estimator in the class of linear predictors, \mathcal{M}=\big\lbrace m:~m(\mathbf{x})=s_h(\mathbf{x})^T\mathbf{y} \text{ where }S=(s(\mathbf{x}_1),\cdots,s(\mathbf{x}_n))^T\text{ is some smoothing matrix}\big\rbrace and assume that y=m_0 (x)+\varepsilon, with \mathbb{E}[\varepsilon]=0 and Var[\varepsilon]=\sigma^2\mathbb{I}, so that m_0 (x)=\mathbb{E}[Y|X=x] . From a theoretical point of view, the quadratic risk, associated with an estimated model \widehat{m}, \mathbb{E}\big[(Y-\widehat{m}(\mathbf{X}))^2\big], is written\mathcal{R}(\widehat{m})=\underbrace{\mathbb{E}\big[(Y-m_0(\mathbf{X}))^2\big]}_{\text{error}}+\underbrace{\mathbb{E}\big[(m_0(\mathbf {X})-\mathbb{E}[\widehat{m}(\mathbf{X})])^2\big]}_{\text{bias}^2}+\underbrace{\mathbb{E}\big[(\mathbb{E}[\widehat{m}(\mathbf{X})]-\widehat{m}(\mathbf{X}))^2\big]}_{\text{variance}} if m_0 is the true model. The first term is sometimes called “Bayes error”, and does not depend on the estimator selected, \widehat{m}.

The empirical quadratic risk, associated with a model m, is here: \widehat{\mathcal{R}}_n(m)=\frac{1}{n}\sum_{i=1}^n (y_i-m(\mathbf{x}_i))^2 (by convention). We recognize here the mean square error, “mse”, which will more generally give the “risk” of the model m when using another loss function (as we will discuss later on). It should be noted that:\displaystyle{\mathbb{E}[\widehat{\mathcal{R}}_n(m)]=\frac{1}{n}\|m_0(\mathbf{x})-m(\mathbf{x})\|^2+\frac{1}{n}\mathbb{E}\big(\|{Y}-m_0(\mathbf{X})\|^2\big)} We can show that:n\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]=\mathbb{E}\big(\|Y-\widehat{m}(\mathbf{x})\|^2\big)=\|(\mathbb{I}-\mathbf{S})m_0\|^2+\sigma^2\|\mathbb{I}-\mathbf{S}\|^2so that the (real) risk of \widehat{m} is: {\mathcal{R}}_n(\widehat{m})=\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]+2\frac{\sigma^2}{n}\text{trace}(\boldsymbol{S})So, if \text{trace}(\boldsymbol{S})\geq0 (which is not a too strong assumption), the empirical risk underestimates the true risk of the estimator. Actually, we recognize here the number of degrees of freedom of the model, the right-hand term corresponding to Mallow’s C_p, introduced in Mallows (1973) using not deviance but R^2.

Statistical Tests

The most traditional test in econometrics is probably the significance test, corresponding to the nullity of a coefficient in a linear regression model. Formally, it is the test of H_0:\beta_k=0 against H_1:\beta_k\neq 0. The so-called Student test, based on the statistics t_k=\widehat{\beta}_k/se_{\widehat{β}_k}, allows to decide between the two alternatives, using the test p-value, defined by \mathbb{P}[|T|>|t_k|] avec T\overset{\mathcal{L}}{\sim} Std_\nu, where \nu is the number of degrees of freedom of the model (\nu=p+1 for the standard linear model). In large dimension, however, this statistic is of very limited interest, given a significant FDR (“False Discovery Ratio”). Classically, with a level of significance \alpha=0.05, 5% of the variables are falsely significant. Suppose that we have p=100 explanatory variables, but that 5 (only) are really significant. We can hope that these 5 variables will pass the Student test, but we can also expect that 5 additional variables (false positive test) will emerge. We will then have 10 variables perceived as significant, while only half are significant, i.e. an FDR ratio of 50%. In order to avoid this recurrent pitfall in multiple tests, it is natural to use the procedure of Benjamini & Hochberg (1995).

From a correlation to some causal effect

Econometric models are used to implement public policy evaluations. It is therefore essential to fully understand the underlying mechanisms in order to know which variables actually make it possible to act on a variable of interest. But then we move on to another important dimension of econometrics. Jerry Neyman was responsible for the first work on the identification of causal mechanisms, and then Rubin (1974) formalized the test, called the “Rubin causal model” in Holland (1986). The first approaches to the notion of causality in econometrics were based on the use of instrumental variables, models with discontinuity of regression, analysis of differences in differences, and natural or unnatural experiments. Causality is usually inferred by comparing the effect of a policy – or more generally of a treatment – with its counterfactual, ideally given by a random control group. The causal effect of the treatment is then defined as \Delta=y_1-y_0, i.e. the difference between what the situation would be with treatment (noted t=1) and without treatment (noted t=0). The concern is that only y=t\cdot y_1+(1-t)\cdot y_0 and t are observed. In other words, the causal effect of variable t  on t  is not observed (since only one of the two potential variables – y_0 or y_1  is observed for each individual), but it is also individual, and therefore a function of x-covariates. Generally, by making assumptions about the distribution of the triplet (Y_0,Y_1,T) , some parameters of the causal effect distribution become identifiable, based on the density of the observable variables (Y,T) . Classically, we will be interested in the moments of this distribution, in particular the average effect of treatment in the population, \mathbb{E}[\Delta] , or even just the average effect of treatment in the case of treatment \mathbb{E}[\Delta|T=1] . If the result (Y_0,Y_1) is independent of the processing access variable T, it can be shown that \mathbb{E}[\Delta]=\mathbb{E}[Y|T=1]- \mathbb{E} [Y|T=0]. But if this independence hypothesis is not verified, there is a selection bias, often associated with \mathbb{E}[Y_0|T=1]- \mathbb{E} [Y_0|T=0]. Rosenbaum & Rubin (1983) propose to use a propensity to be treated score, p(x)=\mathbb{P}[T=1|X=x] , noting that if variable Y_0\ is independent of access to treatment T conditionally to the explanatory variables X, then it is independent of T  conditionally to the score p(X) : it is sufficient to match them using their propensity score. Heckman et al (2003) thus proposes a kernel estimator on the propensity score, which simply provides an estimator of the effect of the treatment, provided that it is treated.

To be continued next time, we’ll introduce “machine learning techniques” (references mentioned above are online here)

Régression sur une variable qualitative et ANOVA

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

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=data.frame(Y=Davis$weight * 2.204622,
                 X1=Davis$sex)

On voulait estimer le modèle y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_iOn avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') 
 Y=Davis$Y

car la matrice \mathbf{X}^T\mathbf{X} est inversible (une fois que l’on enlève la constante)

 solve(t(X)%*%X)
            [,1]       [,2]
[1,] 0.008928571 0.00000000
[2,] 0.000000000 0.01136364

et donc l’estimateur par moindres carrés est (classiquement)\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

 solve(t(X)%*%X) %*% (t(X)%*%Y)
         [,1]
[1,] 125.4272
[2,] 167.3258

ce qui correspond effectivement à la sortie de R,

 reg=lm(Y~0+X1,data=Davis)
 summary(reg)
 
Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X1F  125.427      1.960   64.00   <2e-16 ***
X1M  167.326      2.211   75.68   <2e-16 ***

Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes

x=Y[X[,1]==1]
y=Y[X[,2]==1]
nx=length(x)
ny=length(y)

On avait vu en cours que les \widehat{\mathbf{\beta}} avaient une interprétation très simple, puisque\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_iautrement dit \widehat{{\beta}}_M   est le poids moyen des hommes. Et en effet

 mean(y)
[1] 167.3258

C’est finalement très naturel, ou intuitif.

On peut maintenant s’interroger sur l’écart-type de l’estimateur de \widehat{{\beta}}_M . Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici

 sqrt(var(y)/ny)
[1] 2.794391
 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)
[1] 2.794391

car pour rappel\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de \mathbf{\beta} est proportionnel à \sigma^2 , la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de \sigma^2

 s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2))
 sqrt(s2/ny)
[1] 2.210863

et en effet, on retombe sur la valeur donnée dans le tableau de régresion

 sqrt(s2/nx)
[1] 1.959721

(pareil pour l’autre coefficient).

On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).

 X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y)
          [,1]
[1,] 125.42724
[2,]  41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)
[1] 41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  125.427      1.960   64.00   <2e-16 ***
X1M           41.899      2.954   14.18   <2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -6.4475, df = 286, p-value = 4.826e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -30.62603 -16.30035
sample estimates:
mean of x mean of y 
 143.8626  167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable \mathbf{y} ), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Call:
lm(formula = Y ~ X1, data = Davis)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-57.458 -22.184  -5.512  17.809 118.912 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 143.4382     2.7820   51.56   <2e-16 ***
X1M           0.9645     4.1940    0.23    0.818    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 29.44 on 198 degrees of freedom
Multiple R-squared:  0.000267,	Adjusted R-squared:  -0.004782 
F-statistic: 0.05289 on 1 and 198 DF,  p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y
 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -0.22998, df = 198, p-value = 0.8183
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.235209  7.306165
sample estimates:
mean of x mean of y 
 143.4382  144.4027

et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents.

Copulas and Financial Time Series

I was recently asked to write a survey on copulas for financial time series. The paper is, so far, unfortunately, in French, and is available on https://hal.archives-ouvertes.fr/. There is a description of various models, including some graphs and statistical outputs, obtained from read data.

To illustrate, I’ve been using weekly log-returns of (crude) oil prices, Brent, Dubaï and Maya.

The dataset is available from an excel file, oil.xls (I thought it was possible to load it direclty from the internet, but it did not work… so I suggest to download the file first, and then load it)

> library(xlsx)
> temp <- tempfile()
> download.file(
+ "http://freakonometrics.free.fr/oil.xls",temp)
trying URL 'http://freakonometrics.free.fr/oil.xls'
Content type 'application/vnd.ms-excel' length 99328 bytes (97 KB)
downloaded 97 KB
> oil=read.xlsx(temp,sheetName="DATA",dec=",")
Error in .jcall("RJavaTools", "Ljava/lang/Object;", "invokeMethod", cl,  : 
  java.io.IOException: block[ 0 ] already removed - does your POIFS have circular or duplicate block references?
> oil=read.xlsx("D:\\home\\acharpen\\mes documents\\oil.xls",sheetName="DATA")

Then we can plot those three time series

> head(oil)
        Date      WTI    brent   Dubai     Maya
1 1997-01-10  2.73672  2.25465  3.3673   1.5400
2 1997-01-17 -3.40326 -6.01433 -3.8249  -4.1076
3 1997-01-24 -4.09531 -1.43076 -6.6375  -4.6166
4 1997-01-31 -0.65789  0.34873  0.7326  -1.5122
5 1997-02-07 -3.14293 -1.97765 -0.7326  -1.8798
6 1997-02-14 -5.60321 -7.84534 -7.6372 -11.0549

> Time=as.Date(oil$Date,"%Y-%m-%d")
> plot(Time,oil[,3],type="l",ylab="Brent, weekly log returns",ylim=range(oil[,3:5]))

The idea is to use some multivariate ARMA-GARCH processes here. The heuristics here is that the first part is used to model the dynamics of the average value of the time series, and the second part is used to model the dynamics of the variance of the time series. Two kinds of models are considered in the paper

  • a mutivariate GARCH process (or a model on the dynamics of the variance matrix) on the residuals from the ARMA models
  • a multivariate model (based on copulas) on the residuals of the ARMA-GARCH process

Continue reading Copulas and Financial Time Series

Conditional Distributions from some Elliptical Vectors

This winter, in my ACT8595 course, I asked my students (that was some homework) to prove that it was possible to derive the conditional distribution when we have a Student-t random vector (and to get the analytical expression of the later). But before, let us recall a standard result about the Gaussian vector. If  is a Gaussian random vector, i.e.

then  has a Gaussian distribution. More precisely, it is a  distribution, with

and  is the Schur complement of the block  of the matrix ,

Observe that  is also related to well known quantity: in the bivariate case, where  and  are univariate Gaussian variables,

which is the slope in the linear regression of  on .

In the case of the Student-t distribution, the conditional distrubution will not be a Student-t distribution anymore, but it will still be an elliptical distribution, and some interpretations of various quantities can actually be obtained.

The density of the multivariate centred Student-t distribution, with unit variance, and parameters  and  is

https://latex.codecogs.com/gif.latex?f(\boldsymbol{x})=%20\frac{\Gamma([d+\nu]/2)}{(\nu\pi)^{d/2}%20\Gamma(\nu/2)\vert\boldsymbol{R}\vert^{1/2}}%20\left(%201+\frac{1}{\nu}\boldsymbol{x}%27\boldsymbol{R}^{-1}\boldsymbol{x}%20\right)^{-(d+\nu)/2}

If we consider the following blocks,

https://latex.codecogs.com/gif.latex?\boldsymbol{R}=%20\left(%20\begin{array}{cc}%20\boldsymbol{R}_{11}&%20\boldsymbol{R}_{12}\\%20\boldsymbol{R}_{21}&%20\boldsymbol{R}_{22}%20\end{array}%20\right)

then we can get that marginal distributions have a centred Student-t distribution, with unit variance, and parameters  and ,

https://latex.codecogs.com/gif.latex?f_2(\boldsymbol{x}_2)=%20\frac{\Gamma([d_2+\nu]/2)}{(\nu\pi)^{d_2/2}%20\Gamma(\nu/2)\vert\boldsymbol{R}_{22}\vert^{1/2}}%20\left(%201+\frac{1}{\nu}\boldsymbol{x}_2%27\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_2%20\right)^{-(d_2+\nu)/2}

Then, to derive the conditional density, we can use Bayes formula,

https://latex.codecogs.com/gif.latex?f_{1\vert%202}(\boldsymbol{x}_1\vert%20\boldsymbol{x}_2)=%20\frac{f(\boldsymbol{x}_1,\boldsymbol{x}_2)}{f_2(\boldsymbol{x}_2)}

One can write (as in Section 9.1 in Tong, 1990, The Multivariate Normal Distribution, but other expressions can be found in Section 2.5 in Fang, Ng and Kotz, 1989, Symmetric multivariate and related distributions, or in Section 1.11 in Kotz and Nadarajah, 2004, Multivariate t distributions and their applications) this conditional density as

https://latex.codecogs.com/gif.latex?f_{1\vert%202}(\boldsymbol{x}_1\vert%20\boldsymbol{x}_2)=\kappa%20\left(1+\frac{1}{\nu}\boldsymbol{x}_2%27\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_2\right)^{(d_2+\nu)/2}%20\left(1+\frac{1}{\nu}\left[\boldsymbol{x}_2%27\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_2+\alpha(\boldsymbol{x}_1,\boldsymbol{x}_2)\right]\right)^{-(d_1+\nu)/2}

with

https://latex.codecogs.com/gif.latex?\kappa=\frac{\Gamma([d+\nu]/2)}{(\nu\pi)^{d_1/2}%20\Gamma([d_2+\nu]/2)}\frac{1}{\vert\boldsymbol{R}_{11}-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{R}_{21}\vert^{1/2}}

and

https://latex.codecogs.com/gif.latex?\alpha(\boldsymbol{x}_1,\boldsymbol{x}_2)=(\boldsymbol{x}_1-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_{2})%27%20[\boldsymbol{R}_{11}-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{R}_{21}]^{-1}(\boldsymbol{x}_1-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_{2})

This conditional distribution is elliptical, but it is not a Student-t distribution, except in the case where , or when the correlation matrix  is the identity.

[June 2016] actually, as shown in Ding (2016), this is a Student-t distribution. “Kotz & Nadarajah (2004) and Nadarajah & Kotz (2005) failed to recognize that the conditional distribution of the MVT distribution is also a MVT distribution due to the complexity of the conditional density function […] Conditional distributions of elliptically contoured distributions are also elliptically contoured distributions. But this does not immediately guarantee that conditional distributions of the MVT distributions are also MVT distributions without some further algebra.

Now, if we look at the components of this density, we can observe that we have

https://latex.codecogs.com/gif.latex?(\boldsymbol{x}_1-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_{2})

which was mentioned previously, in the Gaussian case: the term on the right is the conditional mean,

and the bloc that appears at several places is the conditional variance,

Now, if we want to visualize that conditional density, let us plot it. The code below is based on Bayes formula

> library(mnormt)
> r=.6
> R=matrix(c(1,r,r,1),2,2)
> nu=4
> f2=function(x2) dt(x2,df=nu)
> f =function(x) dmt(x,S=R,df=nu)
> f1.2=function(x1,x2) f(c(x1,x2))/f2(x2)

In order to compare that conditional density with a Student-t one, let us define the density of a non-centred Student-t random variable,

> dstd=function(x,mu,s,nu) gamma((nu+1)/2)/
+ (gamma(nu/2)*s*sqrt(pi*nu))*
+ (1+1/nu*(x-mu)^2/(s^2))^(-(nu+1)/2)

Here is the function we can use to plot those two densities,

> graphdensity=function(x2=-1.5){
+ vectx1=seq(-3,3,length=251)
+ y=Vectorize(function(x) f1.2(x,x2))(vectx1)
+ plot(vectx1,y,type="l",col="red",ylim=c(0,.5),
+ xlab="",ylab="")
+ abline(v=r*x2,lty=2)
+ lines(vectx1,dstd(vectx1,x2*r,sqrt(1-r^2),nu),col="blue",lty=2)}
> graphdensity(-1.5)

In the case where , the two lines are rather close (the difference migth come from computational issues)

> graphdensity(-1)

and just to conclude, a last one

> graphdensity(0)

Bivariate Densities with N(0,1) Margins

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are  distributions.

  • The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins , and with correlation . This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
contour(vx,vy,z,col="blue",add=TRUE)
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green")
lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
dnorm(density(X[,1])$x),col="red")
par(mar=c(3,0,1,1))
barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0, 
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks[1])/
diff(yhist$breaks)[1],col="red")

That was the simple part.

  • The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair , we would like to consider the pair , so that the marginal distributions are . The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z)

Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

  • The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where  and  with independent random variables, given , and that both variables are exponentially distributed, with parameter . As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be . The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that  is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value))
Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy[2],z)*dexp(xy[1],z)*
dgamma(z,a,b),lower=0,upper=Inf)$value
f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y))
vx=seq(-3,3,length=151)
vy=seq(-3,3,length=151)
z=matrix(NA,length(vx),length(vy))
for(i in 1:length(vx)){
for(j in 1:length(vy)){
z[i,j]=f(vx[i],vy[j])}}
set.seed(1)
Theta=rgamma(1500,a,b)
Z=cbind(rexp(1500,Theta),rexp(1500,Theta))
X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2]))

There is a small technical problem, but no big deal.

Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian.

Happy St Patrick’s Day

I love Saint Patrick’s Day for, at least, two reasons. The first one is that, on March 17th, you can play out loud The Pogues, the second one is that it’s the only day in the year when I really enjoy getting a Guiness in a pub. And Guiness is important in statistical science (I did mention a couple of hours ago – on this blog –  that beers were important for social reasons in the academic world, but that was for other reasons…)

> theta=seq(0,pi/2,length=101)
> leaf=sin(2*theta)+.25*sin(6*theta)
> for(k in 0:3)
+ polygon(leaf*cos(theta+k*pi/2),leaf*sin(theta+k*pi/2),col="green")

As mentioned in all my statistics and econometrics courses, the history of statistics (I mean here mathematical statistics) is closely related to Guinness.

A long time ago, there was a Guinness Brewing Company of Dublin, which – as its name suggests – was an Irish brewing company. And the boss, who was to inherit the family business, decided to attract young students, trained in chemistry at Cambridge or Oxford.

In 1899, William Sealy Gosset, who had obtained a double degree in math and chemistry, left Oxford to Dublin. And to be quite honest, being graduate in maths meant when he had studied differential equations and astronomy. Basically, mathematics were useless for Guinness, and he got there with his expertise in chemistry. In fact, William turned out to be also a very good administrator, but this has nothing to do with our story.

William had good memories of his studies in math, and he wondered if he could find a problem to look at. He started studies on workmanship, noting that conditions vary so much (temperature, from hops, malt, manufacturing conditions …) that there were only few consistent data. The “law of errors”  (the central limit theorem) can not apply under these conditions.

In short, Bill (now we know each other a little, we’ll call him Bill) took many measurements, and noticed that the Poisson distribution could be an interesting model to work with. To make the story short, Bill managed to use statistical techniques to control the variance of the production, meaning that he was able to lower losses in the production of beer.

A nice application like this one deserved publication in a scientific journal … Well, of course the Poisson distribution has long been known (it was 1904 and a few months before, Von Bortkiewicz found elegant applications of this law, as discussed in a post  a few weeks ago). But there was a disclosure issue there: Bill’s contract prohibited him from disclosing secrets to the competitors.

Meanwhile, Bill had met Karl Pearson, who was then editor of Biometrika, and encouraged him to publish his results. In 1906, Bill who had helped Guiness to gain a lot of money – doing applied mathematics can be usefull – managed to take a sabbatical to work with Pearson to Galton Laboratory biometrics. Bill and Karl decided to publish the work under a pseudonym “Student.” The legend claims that they had hesitated to use “pupil.”

And for almost 30 years, “Mr Gosset” honorable employee Guinness led a dissolute life by publishing in statistical journals (after work in the brewery) always under the pseudonym “Student”. Of course, it might not be that simple. I mean, Bill had a family life, too. And his wife was the captain of the national Hockey team. So I hardly imagine Bill playing the smart ass and doing mathematical computations, when it was time to wash the dishes or iron his shirt…

In 1908, he wrote a remarkable “the probable error of the mean” remarked, at least, by Ronald Fisher. In fact, Bill found that there was a interesting law, but – as the normal – it was difficult to manipulate to obtain confidence intervals. Without a computer, he had the idea of ​​using monte carlo methods to tabulate quantiles and construct its tables. And he was probably the first one to look carefully at the problem of small samples, unlike Karl Pearson, who always put focus on the asymptotic case.

In fact, looking at his small sample, he saw the denominator magnitudes very close to those specifically manipulated Karl, in particular a square root of chi-square law. Well, of course, remained the normality assumption, but at least we had some results for finite samples !

For the story, William Gosset suggested to use letter z for its statistics, the ratio between the mean and (empirical) standard deviation. But a few years later, statisticians became accustomed to use this letter for Gaussian distribution (i.e. when the variance is known), and it became the standard to use the letter t. Hence finally the present name of “Student-t distribution” and in regression outputs, we have the “t-test”.

A legend (told by Harold Hotelling in his memoirs) claims that the Guinness family discovered this double life on the day of the death of William Gosset in 1937 when mathematicians requested financial assistance to print a volume of the works of their employee. But another legend claims that Mr Guinness himself would have suggested his nickname when he had expressed his intention to publish his research… So I guess we’ll never know. But at least, I’ll think about Bill when I’ll get my first Guiness tonight (but I will probably not be able to tell this story anymore when I’ll reach the fourth…)

Test, valeur critique et p-value

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

  • Les erreurs dans un test statistique

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

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

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

  • La p-value

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

  • Valeur critique versus p-value

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

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

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

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

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

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

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

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

image manquante

La statistique de test est ici

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

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

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

Numériquement, on a ici

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

et pour les degrés de liberté

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

La valeur critique est obtenue en lisant dans les tables,

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

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

image manquante

On peut aussi calculer la p-value,

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

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

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

Welch Two Sample t-test

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

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

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

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

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

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

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

One-sample Kolmogorov-Smirnov test

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

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

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

image manquante


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

(nonparametric) copula density estimation

Today, we will go further on the inference of copula functions. Some codes (and references) can be found on a previous post, on nonparametric estimators of copula densities (among other related things).  Consider (as before) the loss-ALAE dataset (since we’ve been working a lot on that dataset)

> library(MASS)
> library(evd)
> X=lossalae
> U=cbind(rank(X[,1])/(nrow(X)+1),rank(X[,2])/(nrow(X)+1))

The standard tool to plot nonparametric estimators of densities is to use multivariate kernels. We can look at the density using

> mat1=kde2d(U[,1],U[,2],n=35)
> persp(mat1$x,mat1$y,mat1$z,col="green",
+ shade=TRUE,theta=s*5,
+ xlab="",ylab="",zlab="",zlim=c(0,7))

or level curves (isodensity curves) with more detailed estimators (on grids with shorter steps)

> mat1=kde2d(U[,1],U[,2],n=101)
> image(mat1$x,mat1$y,mat1$z,col=
+ rev(heat.colors(100)),xlab="",ylab="")
> contour(mat1$x,mat1$y,mat1$z,add=
+ TRUE,levels = pretty(c(0,4), 11))

http://freakonometrics.blog.free.fr/public/perso6/3dcop-est1.gif

Kernels are nice, but we clearly observe some border bias, extremely strong in corners (the estimator is 1/4th of what it should be, see another post for more details). Instead of working on sample https://latex.codecogs.com/gif.latex?(U_i,V_i) on the unit square, consider some transformed sample https://latex.codecogs.com/gif.latex?(Q(U_i),Q(V_i)), where https://latex.codecogs.com/gif.latex?Q:(0,1)\rightarrow\mathbb{R} is a given function. E.g. a quantile function of an unbounded distribution, for instance the quantile function of the https://latex.codecogs.com/gif.latex?\mathcal{N}(0,1) distribution. Then, we can estimate the density of the transformed sample, and using the inversion technique, derive an estimator of the density of the initial sample. Since the inverse of a (general) function is not that simple to compute, the code might be a bit slow. But it does work,

> gaussian.kernel.copula.surface <- function (u,v,n) {
+   s=seq(1/(n+1), length=n, by=1/(n+1))
+   mat=matrix(NA,nrow = n, ncol = n)
+ sur=kde2d(qnorm(u),qnorm(v),n=1000,
+ lims = c(-4, 4, -4, 4))
+ su<-sur$z
+ for (i in 1:n) {
+     for (j in 1:n) {
+ 	Xi<-round((qnorm(s[i])+4)*1000/8)+1;
+ 	Yj<-round((qnorm(s[j])+4)*1000/8)+1
+ 	mat[i,j]<-su[Xi,Yj]/(dnorm(qnorm(s[i]))*
+ 	dnorm(qnorm(s[j])))
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ }

Here, we get

http://freakonometrics.blog.free.fr/public/perso6/3dcop-est2.gif

Note that it is possible to consider another transformation, e.g. the quantile function of a Student-t distribution.

> student.kernel.copula.surface =
+  function (u,v,n,d=4) {
+  s <- seq(1/(n+1), length=n, by=1/(n+1))
+  mat <- matrix(NA,nrow = n, ncol = n)
+ sur<-kde2d(qt(u,df=d),qt(v,df=d),n=5000,
+ lims = c(-8, 8, -8, 8))
+ su<-sur$z
+ for (i in 1:n) {
+     for (j in 1:n) {
+ 	Xi<-round((qt(s[i],df=d)+8)*5000/16)+1;
+ 	Yj<-round((qt(s[j],df=d)+8)*5000/16)+1
+ 	mat[i,j]<-su[Xi,Yj]/(dt(qt(s[i],df=d),df=d)*
+ 	dt(qt(s[j],df=d),df=d))
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ }

Another strategy is to consider kernel that have precisely the unit interval as support. The idea is here to consider the product of Beta kernels, where parameters depend on the location

> beta.kernel.copula.surface=
+  function (u,v,bx=.025,by=.025,n) {
+  s <- seq(1/(n+1), length=n, by=1/(n+1))
+  mat <- matrix(0,nrow = n, ncol = n)
+ for (i in 1:n) {
+     a <- s[i]
+     for (j in 1:n) {
+     b <- s[j]
+ 	mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) *
+     dbeta(b,v/by,(1-v)/by)) / length(u)
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ }

http://freakonometrics.blog.free.fr/public/perso6/3dcop-est3.gif

On those two graphs, we can clearly observe strong tail dependence in the upper (right) corner, that cannot be intuited using a standard kernel estimator…

Copulas estimation and influence of margins

Just a short post to get back on results mentioned at the end of the course. Since copulas are obtained using (univariate) quantile functions in the joint cumulative distribution function, they are – somehow – related to the marginal distribution fitted. In order to illustrate this point, consider an i.i.d. sample http://freakonometrics.blog.free.fr/public/perso6/cop-marg-01.gif from a Student-t distribution,

library(mnormt)
r=.5
n=200
X=rmt(n,mean=c(0,0),S=matrix(c(1,r,r,1),2,2),df=4)

Thus, the true copula is Student-t. Here, with 4 degrees of freedom. Note that we can easily get the (true) value of the copula, on the diagonal

dg=function(t) pmt(qt(t,df=4),mean=c(0,0),
S=matrix(c(1,r,r,1),2,2),df=4)
DG=Vectorize(dg)

Four strategies are considered here to define pseudo-copula base variates,

  • misfit: consider an invalid marginal estimation: we have assumed that margins were Gaussian, i.e. http://freakonometrics.blog.free.fr/public/perso6/cop-marg-2.gif
  • perfect fit: here, we know that margins were Student-t, with 4 degrees of freedom http://freakonometrics.blog.free.fr/public/perso6/cop-marg-3.gif
  • standard fit: then, consider the case where we fit marginal distribution, but in the good family this time (e.g. among Student-t distributions), http://freakonometrics.blog.free.fr/public/perso6/cop-marg-4.gif
  • ranks: finally, we consider nonparametric estimators for marginal distributions, http://freakonometrics.blog.free.fr/public/perso6/cop-marg-10.gif

Now that we have a sample with margins in the unit square, let us construct the empirical copula,

http://freakonometrics.blog.free.fr/public/perso6/cop-marg-6.gif
Let us now compare those four approaches.

  • The first one is to illustrate model error, i.e. what’s going on if we fit distributions, but not in the proper family of parametric distributions.
X0=cbind((X[,1]-mean(X[,1])/sd(X[,1])),
(X[,2]-mean(X[,2])/sd(X[,2])))
Y=pnorm(X0)

Then, the following code is used to compute the value of the empirical copula, on the diagonal,

diagonale=function(t,Z) mean((Z[,1]<=t)&(Z[,2]<=t))
diagY=function(t) diagonale(t,Y)
DiagY=Vectorize(diagY)
u=seq(0,1,by=.005)
dY=DiagY(u)

On the graph below, 1,000 samples of size 200 have been generated. All trajectories are the estimation of the copula on the diagonal. The black plain line is the true value of the copula

Obviously, it is not good at all. Mainly because the distribution of http://freakonometrics.blog.free.fr/public/perso6/cop-marg-8.gif can’t be a copula, since margins are not even uniform on the unit interval.

  • a perfect fit. Here, we use the following code to generate our copula-type sample
U=pt(X,df=4)

This time, the fit is much better.

  • Using maximum likelihood estimators to fit the best distribution within the Student-t family
F1=fitdistr(X0[,1],dt,list(df=5),lower = 0.001)
F2=fitdistr(X0[,2],dt,list(df=5),lower = 0.001)
V=cbind(pt(X0[,1],df=F1$estimate),pt(X0[,2],df=F2$estimate))

Here, it is also very good. Even better than before, when the true distribution is considered.

(it is like using Lillie test for goodness of fit, versus Kolmogorov-Smirnov, seehere for instance, in French).

  • Finally, let us consider ranks, or nonparametric estimators for marginal distributions,
R=cbind(rank(X[,1])/(n+1),rank(X[,2])/(n+1))

Here it is even better then the previous one

If we compare Box-plots of the value of the copula at point (.2,.2), we obtain the following, with on top ranks, then fitting with the good family, then using the true distribution, and finally, using a non-proper distribution.

Just to illustrate one more time a result mentioned in a previous post, “in statistics, having too much information might not be a good thing“.

Does the Student based confidence interval have any interest in practice ?

Friday in the course of statistics, we started the section on confidence interval, and like always, I got a bit confused with the degrees of freedom of the Student (should it be http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif or http://freakonometrics.blog.free.fr/public/perso2/IC-std-5.gif ?) and which empirical variance (should we consider the one where we divide by http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif or the one with http://freakonometrics.blog.free.fr/public/perso2/IC-std-5.gif ?).
And each time I start to get confused, the student obviously see it, and start to ask tricky questions… So let us make it clear now. The correct formula is the following: let

http://freakonometrics.blog.free.fr/public/perso2/IC-std-4.gif

then

http://freakonometrics.blog.free.fr/public/perso2/IC-std-1.gif

is a confidence interval for the mean of a Gaussian i.i.d. sample.
But the important thing is neither the n-1 that appear as degrees of freedom nor the http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif that appear in the estimation of the standard error. Like always in mathematical result, the most important part of that result is not mentioned here: observations have to be i.i.d. and to be normally distributed. And not “almost” normally distributed….
Consider the following case: we have http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif=20 observations that are almost normally distributed. Hence, I consider a student t distribution

n=20; X=rt(n,df=3)

An Anderson Darling normality test accepts a normal distribution in 2 cases out of 3.

for(s in 1:10000){
X=rt(n,df=3)
pv[s]=ad.test(X)$p.value
}
mean(pv>.05)
[1] 0.6799

With a true normal distribution if would be 95% of the cases, so in some sense, I can pretend that I generate almost normal samples.
For those samples, we can look at bounds of the 90% confidence interval for the mean, with three different formulas,

http://freakonometrics.blog.free.fr/public/perso2/IC-std-1.gif

i.e. the correct one, or the one where I considered http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif degrees of freedom instead of http://freakonometrics.blog.free.fr/public/perso2/IC-std-5.gif,

http://freakonometrics.blog.free.fr/public/perso2/IC-std-2.gif

and the one were we condired a Gaussian quantile instead of a Student t one,

http://freakonometrics.blog.free.fr/public/perso2/IC-std-3.gif

(and one might think to look at the non-unbiased estimator of the variance, also).
for(s in 1:10000){
X=rt(n,df=3)
m[s]=mean(X)
sd=sqrt(var(X))
IC1[s]=m[s]-qt(.95,df=n-1)*sd/sqrt(n)
IC2[s]=m[s]-qt(.95,df=n)*sd/sqrt(n)
IC3[s]=m[s]-qnorm(.95)*sd/sqrt(n)
}

One the graph below are plotted the distributions of the values obtained as lower bound of the 90% confidence interval,

(the curves with http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif and http://freakonometrics.blog.free.fr/public/perso2/IC-std-5.gif degrees of freedom in quantiles are the same, here).
The dotted vertical line is the true lower bound of the 90%-confidence interval, given the true distribution (which was not a Gaussian one).
If I get back to the standard procedure in any statistical textbook, since the sample is almost Gaussian, the lower bound of the confidence interval should be (since we have a Student t distribution)

mean(IC1)
[1] -0.605381

instead of

mean(IC3)
[1] -0.5759391

(obtained with a Gaussian distribution instead of a Student one). Actually, both of them are quite different from the correct one which was

quantile(m,.05)
       5% 
-0.623578

As I mentioned in a previous post (here), an important issue is that if we do not know a parameter and substitute an estimator, there is usually a cost (which means usually that the confidence interval should be larger). And this is what we observe here. From a teacher’s point of view, it is an important issue that should be mentioned in statistical courses….

But another important point is also that confidence interval is valid only if the underlying distribution is Gaussian. And not almost Gaussian, but really a Gaussian one.  So since with http://freakonometrics.blog.free.fr/public/perso2/IC-std-6.gif=20 observations everything might look Gaussian, I was wondering what should be done in practice… Because in some sense, using a Student quantile based confidence interval on some almost Gaussian sample is as wrong as using a Gaussian quantile based confidence interval on some Gaussian sample…

In statistics, having too much information might not be a good thing

A common idea in statistics is that if we don’t know something, and we use anestimator of that something (instead of the true value) then there will be some additional uncertainty. For instance, consider a random sample, i.i.d., from a Gaussian distribution. Then, a confidence interval for the mean is

http://freakonometrics.blog.free.fr/public/perso2/IC-cout-06.gif

where http://freakonometrics.blog.free.fr/public/perso2/inc-out-8.gif is the quantile of probability level http://freakonometrics.blog.free.fr/public/perso2/IC-cout-05.gif of the standard normal distribution http://freakonometrics.blog.free.fr/public/perso2/inc-out-09.gif. But usually, standard deviation http://freakonometrics.blog.free.fr/public/perso2/inc-cout-10.gif (the something is was talking about earlier) is usually unknown. So we substitute an estimation of the standard deviation, e.g.

http://freakonometrics.blog.free.fr/public/perso2/IC-cout-02.gif

and the cost we have to pay is that the new confidence interval is

http://freakonometrics.blog.free.fr/public/perso2/IC-cout-01.gif

where now http://freakonometrics.blog.free.fr/public/perso2/IC-cout-03.gif is the quantile of the Student distribution, of probability level http://freakonometrics.blog.free.fr/public/perso2/IC-cout-05.gif, with http://freakonometrics.blog.free.fr/public/perso2/IC-cout-04.gif degrees of freedom.
We call it a cost since the new confidence interval is now larger (the Student distribution has higher upper-quantiles than the Gaussian distribution).
So usually, if we substitute an estimation to the true value, there is a price to pay.
A few years ago, with Jean David Fermanian and Olivier Scaillet, we were writing a survey on copula density estimation (using kernels,  here). At the end, we wanted to add a small paragraph on the fact that we assumed that we wanted to fit a copula on a sample http://freakonometrics.blog.free.fr/public/perso2/ic-cout_11.gif i.i.d. with distribution http://freakonometrics.blog.free.fr/public/perso2/ic-cout_13.gif, a copula, but in practice, we start from a samplehttp://freakonometrics.blog.free.fr/public/perso2/ic-cout_12.gif with joint distribution http://freakonometrics.blog.free.fr/public/perso2/ic-cour_14.gif (assumed to have continuous margins, and – unique – copula http://freakonometrics.blog.free.fr/public/perso2/ic-cout_13.gif). But since margins are usually unknown, there should be a price for not observing them.
To be more formal, in a perfect wold, we would consider

http://freakonometrics.blog.free.fr/public/perso2/ic-cout-15.gif

but in the real world, we have to consider

http://freakonometrics.blog.free.fr/public/perso2/ic-cout-16.gif

where it is standard to consider ranks, i.e. http://freakonometrics.blog.free.fr/public/perso2/ic-cout_109.gif are empirical cumulative distribution functions.
My point is that when I ran simulations for the survey (the idea was more to give illustrations of several techniques of estimation, rather than proofs of technical theorems) we observed that the price to pay… was negative ! I.e. the variance of the estimator of the density (wherever on the unit square) was smaller on the pseudo sample http://freakonometrics.blog.free.fr/public/perso2/ic-cout-17.gif than on perfect sample http://freakonometrics.blog.free.fr/public/perso2/ic-cout_18.gif.
By that time, we could not understand why we got that counter-intuitive result: even if we do know the true distribution, it is better not to use it, and to use instead a nonparametric estimator. Our interpretation was based on the discrepancy concept and was related to the latin hypercube construction:

With ranks, the data are more regular, and marginal distributions are exactlyuniform on the unit interval. So there is less variance.
This was our heuristic interpretation.
A couple of weeks ago, Christian Genest and Johan Segers proved that intuition in an article published in JMVA,

Well, we observed something for finite http://freakonometrics.blog.free.fr/public/maths/mariage01.png, but Christian and Johan obtained an analytical result. Hence, if we denote

http://freakonometrics.blog.free.fr/public/perso2/JSCG-1.gif

the empirical copula in the perfect world (with known margins) and

http://freakonometrics.blog.free.fr/public/perso2/JSCG-2.gif

the one constructed from the pseudo sample, they obtained that, everywhere

http://freakonometrics.blog.free.fr/public/perso2/JSCG-6.gif

with nice graphs of http://freakonometrics.blog.free.fr/public/perso2/JSCG-7.gif,

So I was very happy last week when Christian show me their results, to learn that our intuition was correct. Nevertheless, it is still a very counter-intuitive result…. If anyone has seen similar things, I’d be glad to hear about it !

Happy birthday Arthur !

mais non, ce n’est pas mon anniversaire ! Mais l’occasion était trop belle…

oui, c’est l’anniversaire d’Arthur Guinness (Arth Guinness pour ceux qui ont déjà eu une canette entre les mains…)

Et comme je l’ai rappelé lors des rappels de statistiques, l’histoire de la statistique (j’entends par là statistique mathématique) est étroitement liée à Guiness.

Il y fort fort longtemps, il existait une Guiness Brewing Company of Dublin, qui comme son nom l’indique, était une brasserie irlandaise. Et au début du XXème siècle, le propriétaire, qui venait d’hériter de cette entreprise familiale décida d’attirer des jeunes étudiants, formés en chimie à Cambridge, ou Oxford.

En 1899, un  certain William Sealy Gosset, qui avait obtenu un double diplôme, en maths et en chimie, quitta Oxford pour Dublin. Et pour être tout à fait honnête, diplômé en maths signifiait alors qu’il avait étudié les équations différentielles ou l’astronomie. En gros, les maths, ça ne servait à rien chez Guiness, et il était arrivé là grâce à ses compétences en chimie. Et en fait, William s’avérait être en plus un très bon administrateur, mais on sort un peu de notre sujet.

Mais William avait gardé un bon souvenir de ses études en maths, et il se demandait s’il pouvait trouver un problème sur lequel se pencher…. il commença a travailler des études de qualité de fabrication, notant que les conditions varient tellement (température, provenance du houblon, du malt, conditions de fabrication…) que les données homogènes sont peu nombreuses. La “loi des erreurs” classique ne peut pas s’appliquer dans ces conditions, i.e. le théorème central limite.

Bref, Bill (maintenant qu’on se connaît un peu, on va l’appeler Bill) pris beaucoup de mesures, et observa une loi de Poisson, ce qui était pratique car cette loi n’avait qu’un paramètre. Bref, Bill arriva à contrôler la variance, et donc à faire moins de pertes lors de la production de bière…

Une jolie application comme celle là méritait publication dans une revue scientifique… Bon, bien sûr la loi de Poisson était connue depuis longtemps (on était en 1904 et quelques mois auparavant, Von Bortkiewicz avait trouvé d’élégantes applications à cette loi, en particulier sur les accidents de chevaux dans l’armée Prusse). Et à côté, il y avait un soucis de taille que connaissent bien des actuaires travaillant dans les compagnies d’assurance: le contrat de Bill lui interdisait de divulguer des secrets à la concurrence, qui plus est des secrets capables de leur faire gagner de l’argent !

Entre temps, Bill avait rencontré Karl Pearson, qui était alors éditeur de la revue Biometrika, et qui l’encouragea à publier ses résultats. En 1906, Bill qui avait fait gagné beaucoup d’argent à Guiness – et qui avait d’ailleurs montré l’intérêt de faire des maths appliqués – réussi à prendre une année sabbatique pour aller travailler avec Pearson au laboratoire de biométrie de Galton. Et Bill et Karl décidèrent de faire publier les travaux sous un pseudonyme, “Student”. La légende prétend qu’ils avaient aussi hésité à prendre “pupil“.

Et pendant presque 30 ans, “Mr Gosset“, honorable employé de Guinness menait une vie dissolue en publiant dans des revues de statistiques (après son travail dans la brasserie1), toujours sous le pseudonyme de “Student“.

En 1908, on lui doit un remarquable “the probable error of the mean“, qui remarque en tous les cas Ronald Fisher. En fait, Bill trouva qu’il existait une loi assez remarquable, mais qui – comme la loi normale – était délicate à manipuler pour obtenir des intervalles de confiance. Sans ordinateur, il eu malgré tout l’idée d’utiliser des méthodes de monte carlo pour tabuler les quantiles, et construire ses tables. Et il fut un des premiers à se pencher sur le problème des petits échantillons, contrairement à Karl Pearson, qui se plaçait toujours dans le cas asymptotique.

En fait, en regardant son petit échantillon, il vit apparaître au dénominateur des grandeurs très proches de celles que manipulait précisément Karl, en particulier une racine carré de loi du chi-deux. Bon, bien sûr, restait l’hypothèse (fondamentale) de normalité des observations, mais au moins, on avait des résultats à distance finie !

Pour la petite histoire, William Gosset proposa la lettre z pour sa statistique, le ratio entre la moyenne et l’écart-type. Mais quelques années plus tard, les statisticiens prirent l’habitude d’utiliser cette lettre pour des lois normales (autrement dit quand la variance est connue), et la norme devint d’utiliser la lettre t. D’où finalement la dénomination actuelle de “Student-t distribution”.

Une légende (racontée par Harold Hotelling dans ses mémoires) veut que la famille Guinness découvrit cette double vie le jour de la mort de William Gosset, en 1937, lorsque des mathématiciens demandèrent une aide financière pour imprimer un volume des œuvres de leur employé. Mais une autre légende raconte que c’est Guinness qui lui aurait suggéré son pseudo lorsqu’il aurait fait part de ses intentions de publier ses travaux de recherche…

1 et si si, Bill avait une vie de famille… il avait d’ailleurs une femme qui devait le mener à la baguette, puisqu’elle était à l’époque capitaine de l’équipe nationale de Hockey ! il devait pas faire le malin avec ses calculs de loi quand c’était l’heure de faire la vaisselle….