All posts by Arthur Charpentier

Arthur Charpentier, professor in Montréal, in Actuarial Science. Former professor-assistant at ENSAE Paristech, associate professor at Ecole Polytechnique and assistant professor in Economics at Université de Rennes 1.  Graduated from ENSAE, Master in Mathematical Economics (Paris Dauphine), PhD in Mathematics (KU Leuven), and Fellow of the French Institute of Actuaries.

From multinomial regression to binary classification on some Siamese data

There are two kinds of people in the world: people who think there are two kinds of people in the world and people who don’t

(borrowed from Menand (2018)). Because things are always simpler when we face only binary choice, aren’t they? But consider here the case were multiple options are possible, and let us see if we cannot get back to simpler binary choices. Consider a collection of observations (y_i,\boldsymbol{x}_i) where y_i is some categorical variable, y_i\in\mathcal{A} where \mathcal{A}=\lbrace A_1,\cdots,A_\kappa \rbrace, with \kappa possible categories. Let \mathcal{I}_k=\lbrace i:y_i\in A_k \rbrace.

In a classical multinomial logistic regression, suppose that A_1 is the reference, then \mathbb{P}[Y=A_j|\boldsymbol{X}=\boldsymbol{x}]=\frac{\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_j]}{1+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_2]+\cdots+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]}With a lot a categories, and a small number of observations, inference can be complicated, and non-robust.

  • the Siamese dataset

The name Siamese I use, here, comes from Siamese Networks. Or sort of… As we say in French, it is an « histoire de l’homme qui a vu l’homme qui a vu l’ours » (story of the man who saw the man who saw the bear). A few years ago, a student tried to explain to me the idea of Siamese Networks and this is what I understood. I might be completely wrong, but the idea I got from it did make sense, in my mind at least. That is the story of that blog post…

The idea of the siamese algorithm will be to consider all pairs of observations, (y_i,\boldsymbol{x}_i) and (y_j,\boldsymbol{x}_j) :

  1. \tilde y_{i,j}=\boldsymbol{1}(y_i=y_j) indicating if individuals i and j are in the same category
  2. \tilde{\boldsymbol{x}}_{i,j} is a collection of p-1 variables,
  • \tilde {x}_{k:i,j}={x}_{k:i}-{x}_{k:j} if x_k is continuous, or \tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}| (we can use another metric, e.g. \tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|^2, and this is why I decided to use some GAM model in the logistic regression on the Siamese dataset)
  • \tilde {x}_{k:i,j}=({x}_{k:i},{x}_{k:j})\in\mathcal{X}_k\times\mathcal{X}_k if x_k is a categorical variables (taking values in the set \mathcal{X}_k), or \tilde {x}_{k:i,j}=\boldsymbol{1}({x}_{k:i}\neq{x}_{k:j})\in\{0,1\}

The original dataset was a n\times p matrix, and (if there are no categorical variable), it becomes a n(n-1)/2\times p matrix. The key point is that if the original variable y_i was multinomial, y_{i,j} is now binomial. For instance, if our initial dataset was the following, with two covariates, one continuous and one categorical

its siamese counterpart is the following

  • Classification step

On the dataset (\tilde y_{i,j},\tilde {x}_{k:i,j})_{i,j}, fit a logistic regression, \mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}]=\frac{\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}(or any classification model – CART, random forest, etc). But that is the easy part (unless n is large, because the siamese dataset has (roughly) n^2/2 rows). The difficult task is the prediction

  • Prediction step

Consider a new input variable \boldsymbol{x}_{\cdot}, and define its siamese version, \tilde{\boldsymbol{x}}_{\cdot}=(\tilde{\boldsymbol{x}}_{\cdot,j})_j, i.e. a database with n rows. Then compute
p_{\cdot,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{\cdot,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]} where p_{\cdot,j} is the probability that (y_j,\boldsymbol{x}_{j}) and (y_{\cdot},\boldsymbol{x}_{\cdot}) are in the same category, as well as p_{i,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{i,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}
Let \boldsymbol{p}_{\cdot}=(p_{\cdot,j}), and similarly \boldsymbol{p}_{i}=(p_{i,j}). Then several techniques can be used to predict y_{\cdot}.

  1. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{j=1,\cdots,n}{\text{argmax}}\{p_{\cdot,j}\}: the predicted class is the one of the observation the most likely to be other same class
  2. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{\ell=1,\cdots,k}{\text{argmax}}\{\overline{p}_{\ell}\}, where\overline{p}_{\ell} = \frac{1}{n_{\ell}}\sum_{j\in\mathcal{I}_s} \boldsymbol{1}(y_j =y_{\ell}),\text{ where }\mathcal{I}_s=\lbrace i:p_{\cdot,i}>s\rbraceconsider only probabilities sufficiently high to be considered, and predict the most important class (majority rule)
  3. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{\theta_i\} where \theta_i=\cos(\boldsymbol{p}_{\cdot},\boldsymbol{p}_{i})=\displaystyle{\frac{\boldsymbol{p}_{\cdot}\cdot\boldsymbol{p}_{i}}{\|\boldsymbol{p}_{\cdot}\|\|\boldsymbol{p}_{i}\|}}
  4. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{KL_{\cdot|i}\} and KL_{\cdot|i}=\displaystyle{\sum_{j=1}^n p_{\cdot,j}\log\frac{p_{\cdot,j}}{p_{i,j}}} (but one can select another metric)
  5. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{j\in\mathcal{J}}{\text{argmax}}\{p_{\cdot,j}\} and \mathcal{J} is a sample of k observations, chosen randomly, one in each group (one-shot procedure): the predicted class is the one of the observation the most likely to be other same class

Heuristically, it can be related to some k nearest neighbors strategy: we give the attribute that most neighbors have. The total distance is a weighted sum of the componentwise distances (for the logistic regression).

  • Simulation study

In order to test that technique, let us generate some multinomial model where y has 10 possible labels, with 6 (independent) covariates x_1,\cdots,x_6, and \mathbb{P}[Y=A_k|\boldsymbol{X}=\boldsymbol{x}]\propto \exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k] (where coefficients \boldsymbol{\beta}_k where generated randomly) for k\in\{1,2,\cdots,10\} (there were 10 categories) and with n=700 observations.

n=700
X1=rnorm(n)
X2=rnorm(n)
X3=rnorm(n)
X4=rnorm(n)
X5=rnorm(n)
X6=rnorm(n)
X=cbind(1,X1,X2,X3,sqrt(abs(X4)),X5*X1,X6)
k = 10
 PARAM = matrix(rnorm(k*6),k,6)
 PARAM[,1]=PARAM[,1]-1
 PARAM=cbind(PARAM,0)
 P=matrix(NA,n,k-1)
 for(j in 1:(k-1)) P[,j] = X %*% (PARAM[j,])+rnorm(n)
 P=cbind(P,0)
S=apply(exp(P),1,sum)
Pb = exp(P)/S
tirage = function(i){
      sample(1:10,size=1,prob = Pb[i,])
}
Y = LETTERS[Vectorize(tirage)(1:n)]
dbase = data.frame(Y=as.factor(Y),X1,X2,X3,X4,X5)

In the paragraph previously, I suggested to take the most likely one. Being wrong means that it was not the first choice. But perhaps being the second or the third choice is not that bad, actually. So in my simulations, I look at the proportion of predictions were our prediction is the good one (top 1), if the true one is either the most likely or the second most likely (top 2), or in the top 3. That will be on my x-axis. I draw some line, but we simply have three points (top 1, top 2 and top 3). I compute the proportion of good prediction, using cross-validation techniques (10-fold). The black lines are one of the methods described above. The red one is the standard multinomial model (with a logistic link function). For the Siamese model, I tried several models. I tried a logistic regression, and some smooth version (GAM) on top

and a classification tree, on the left, as well as some random forest on the right, below.

It looks like the multinomial approach performs always better than any Siamese one… and to be honest, I am disappointed.

Here is the code I did use when I considered a logistic regression on the Siamese dataset,

set.seed(1)
kfold = sample(rep(1:10,n/10))
 
KFOLDglm = function(i){
i_test=which(kfold==i)
i_calibration=which(kfold!=i)
y=credit[i_calibration,"Y"]
tirage = function(){
v=c(sample(i_calibration[y==levels(y)[1]],size=1),
    sample(i_calibration[y==levels(y)[2]],size=1),
    sample(i_calibration[y==levels(y)[3]],size=1),
    sample(i_calibration[y==levels(y)[4]],size=1),
    sample(i_calibration[y==levels(y)[5]],size=1),
    sample(i_calibration[y==levels(y)[6]],size=1),
    sample(i_calibration[y==levels(y)[7]],size=1),
    sample(i_calibration[y==levels(y)[8]],size=1),
    sample(i_calibration[y==levels(y)[9]],size=1),
    sample(i_calibration[y==levels(y)[10]],size=1))
names(v)=levels(y)
return(v)
}
 
LogisticModel <- multinom(Y ~ ., data = credit[i_calibration,], trace=FALSE)
 
comparaisonx = function(base,x=base[1,]){
  mix_base = base
  for(j in 1:ncol(base)){
    xj = as.numeric(x[j])-base[,j]
    mix_base[,j] = (xj)
  }
  mix_base
}
comparaisony = function(base,y=base[1]){
  as.factor(base == y)
}
creditx = credit[,-which(names(credit) == "Y")]
nc=length(i_calibration)
B=comparaisonx(base = creditx[i_calibration[2:nc],],x=creditx[i_calibration[1],])
B$Y=comparaisony(base = credit[i_calibration[2:nc],"Y"],y=credit[i_calibration[1],"Y"])
for(i in 2:(nc-1)){
  B0=comparaisonx(base = creditx[i_calibration[(i+1):nc],],x=creditx[i_calibration[i],])
  B0$Y=comparaisony(base = credit[i_calibration[(i+1):nc],"Y"],y=credit[i_calibration[i],"Y"])
  B=rbind(B,B0)
}
credit_mix = B
 
OneShotLogisticModel <- glm(Y ~ ., data = credit_mix, family=binomial)
A_ref = table(credit[i_calibration,"Y"])/length(i_calibration)
 
vect_oneshot = function(i){
  B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,])
  predict(OneShotLogisticModel,type="response",newdata=B2)
}
 
prediction_oneshot = function(i,type=1){
B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,])
p=predict(OneShotLogisticModel,type="response",newdata=B2)
y=credit[i_calibration,"Y"]
base = data.frame(p,y)
base = base[rev(order(base$p)),]
if(type==1){T = table(base$y[1:11])
return(names(which.max(T)))}
if(type==2){return(base$y[1])}
if(type==3){A=table(base$y[1:10])/10
T=A/A_ref
return(names(which.max(T)))}
if(type==4){
  costheta = rep(NA,length(i_calibration))
  for(j in 1:length(i_calibration)){
    vecteur_proba = vect_oneshot(i_calibration[j])
    costheta[j] = sum(vecteur_proba*p)/(sqrt(sum(vecteur_proba^2))*sqrt(sum(p^2)))
  }
  return(y[which.max(costheta)])}
if(type==5){
  kl = rep(NA,length(i_calibration))
  for(j in 1:length(i_calibration)){
    vecteur_proba = vect_oneshot(i_calibration[j])
    kl[j] = sum(p*log(vecteur_proba/p))
  }
  return(y[which.max(as.vector(kl))])}
if(type==6){ ## one shot : tirer au hasard un de chaque, et dire lequel est plus credible !
 
y=credit[i_calibration,"Y"]
tirage = function(){
    v=c(sample(i_calibration[y==levels(y)[1]],size=1),
        sample(i_calibration[y==levels(y)[2]],size=1),
        sample(i_calibration[y==levels(y)[3]],size=1),
        sample(i_calibration[y==levels(y)[4]],size=1),
        sample(i_calibration[y==levels(y)[5]],size=1),
        sample(i_calibration[y==levels(y)[6]],size=1),
        sample(i_calibration[y==levels(y)[7]],size=1),
        sample(i_calibration[y==levels(y)[8]],size=1),
        sample(i_calibration[y==levels(y)[9]],size=1),
        sample(i_calibration[y==levels(y)[10]],size=1))
    names(v)=levels(y)
    return(v)
  }
  pd=rep(NA,101)
for(ix in 1:101){
ids = tirage()
B2=comparaisonx(base = creditx[ids,],x=creditx[i,])
p=predict(OneShotLogisticModel,type="response",newdata=B2)
pd[ix]=levels(y)[which.max(p)]
}
levels(y)[which.max(table(pd))]
}
}
PRED0=as.character(credit[i_test,"Y"])
PRED1=as.character(predict(LogisticModel,type = "class",
                           newdata=credit[i_test,]))
PRED21=as.character(Vectorize(function(i) prediction_oneshot(i,type=1))
                    (i_test))
PRED22=as.character(Vectorize(function(i) prediction_oneshot(i,type=2))(i_test))
PRED23=as.character(Vectorize(function(i) prediction_oneshot(i,type=3))(i_test))
PRED24=as.character(Vectorize(function(i) prediction_oneshot(i,type=4))(i_test))
PRED25=as.character(Vectorize(function(i) prediction_oneshot(i,type=5))(i_test))
PRED26=as.character(Vectorize(function(i) prediction_oneshot(i,type=6))(i_test))
B=data.frame(PRED0,PRED1,PRED21,PRED22,PRED23,PRED24,PRED25,PRED26)
B
}
 
s=1/100;setTxtProgressBar(pb, s*2)
for(i in 2:10){
  PREDICTION = rbind(PREDICTION,KFOLDglm(i))
s=s+1/100;setTxtProgressBar(pb, s*2)}
for(j in 1:5) PREDICTION[,j]=as.character(PREDICTION[,j])
L=list()
v=mean(PREDICTION[,1]!=PREDICTION[,2])
names(v)="logistic"
L[["logistic"]]=v
v=c(mean(PREDICTION[,1]!=PREDICTION[,3]),
  mean(PREDICTION[,1]!=PREDICTION[,4]),
  mean(PREDICTION[,1]!=PREDICTION[,5]),
  mean(PREDICTION[,1]!=PREDICTION[,6]),
  mean(PREDICTION[,1]!=PREDICTION[,7]),
  mean(PREDICTION[,1]!=PREDICTION[,8]))
names(v)=c("top10","max","10norm","cos","KL","OS")
L[["glm"]]=v

Autocalibration for Insurance Pricing with Machine Learning

With Michel Denuit and Julien Trufin, we recently uploaded a joint paper on ArXiv, entitled Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning.

Boosting techniques and neural networks are particularly effective machine learning methods for insurance pricing.
Often in practice, there are nevertheless endless debates about the choice of the right loss function to be used to train the machine learning model, as well as about the appropriate metric to assess the performances of competing models. Also,
the sum of fitted values can depart from the observed totals to a large extent and this often confuses actuarial analysts.
The lack of balance inherent to training models by minimizing deviance outside the familiar GLM with canonical
link setting has been empirically documented in Wüthrich (2019, 2020) who attributes it to the early stopping rule in gradient descent methods for model fitting.
The present paper aims to further study this phenomenon when learning proceeds by minimizing Tweedie deviance.
It is shown that minimizing deviance involves a trade-off between the integral of weighted differences of lower partial moments and the bias measured on a specific scale.
Autocalibration is then proposed as a remedy. This new method to correct for bias adds an extra local GLM step to the analysis.
Theoretically, it is shown that it implements the autocalibration concept in pure premium calculation and ensures that balance also holds on a local scale, not only at portfolio level as with existing bias-correction techniques. The convex order appears to be the natural tool to compare competing models, putting a new light on the diagnostic graphs and associated metrics proposed by Denuit et al. (2019).

In this paper, we started with a simple observation : with GLM (like a Poisson regression) we have unbiased predictions, in the sense that \widehat{y}_1+...+\widehat{y}_n=y_1+...+y_n (on the training dataset, we should expect to have \widehat{y}_1+...+\widehat{y}_n\approx y_1+...+y_n on the validation dataset). But with any machin learning algorithm, this is no longer the case. For instance, with a boosting algorithm, we can end up with a significant bias. This is an application on an insurance dataset from the CASdataset package, with a GLM (Poisson regression) on the left, and additive smooth version, GAM in the middle and some boosting algorithm on the right. Here is the dispersion of the predictions, \widehat{y}_i‘s,

Let \widehat{\pi} denote a fited model. If \widehat{\pi} was close to the true parameter of the Poisson model, \mu where \mu(\boldsymbol{x})=\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}], then function  s\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=s] should be very close to the first diagonal (since \mathbb{E}[Y|\mu(\boldsymbol{X})=s]=s). And this is was we observe below, for the GLM. But the boosting algorithm has a significant bias.

Furthermore, the bias is not the same everywhere. We can see it more precisely on a quantile version of the x-axis

that is u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)]. Here is a multiplicative version of that bias u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)]/u

The idea of autocalibration is to use that multiplicative factor to correct from the bias. Thus, if \widehat{\pi}(\boldsymbol{x}) is close to the 75% upper quantile, then the true value should be 20% larger for the bootsting algorithm. Here is the new distribution of the predictions, with that correction

and we can observe that u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)] is now much closer to the first diagonal (below is the multiplicative bias)

We discuss in the paper this correction for local bias, so that u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)] becomes as close as possible to the first diagonal… see https://arxiv.org/abs/2103.03635 for a complete version of the paper, and https://github.com/freakonometrics/autocalibration/ for the complete version of the R codes

Unusual data for insurance, joint research initiative

We recently started a joint research initiative, funded by the AXA Research Fund, to work on unusual data for insurance

Insurers sometimes lack information at the time of claim submission, such as the structure of a building, the presence of health risks common to a group of people, or the spatial diffusion of a pandemic. Unusual data, such as satellite images, personal network connections, and tweets can be used to populate this information gap. In this joint research initiative, we will use images, network data, and texts for risk analysis from an actuarial perspective. Specifically, the project will explore how using unusual data can contribute to smoother claims assessment and reduced data quality risk, allowing for better risk selection and pricing. The project will look at three types of unusual data: pictures/satellite images, network data, and text data.

More information will be shared via a dedicated website (https://jridata.github.io/), even if I will also mention interesting papers, conferences and open-source codes on this blog….

Principal Component Analysis: A Generalized Gini Approach

Our paper, with Stéphane Mussard and Téa Ouraga, entitle Principal Component Analysis: A Generalized Gini Approach is finally out in the European Journal of Operations Research.

A principal component analysis based on the generalized Gini correlation index is proposed (Gini PCA). The Gini PCA generalizes the standard PCA based on the variance. It is shown, in the Gaussian case, that the standard PCA is equivalent to the Gini PCA. It is also proven that the dimensionality reduction based on the generalized Gini correlation matrix, that relies on city-block distances, is robust to outliers. Monte Carlo simulations and an application on cars data (with outliers) show the robustness of the Gini PCA and provide different interpretations of the results compared with the variance PCA.

Some general thoughts on Partial Dependence Plots with correlated covariates

The partial dependence plot is a nice tool to analyse the impact of some explanatory variables when using nonlinear models, such as a random forest, or some gradient boosting.The idea (in dimension 2), given a model m(x_1,x_2) for \mathbb{E}[Y|X_1=x_1,X_2=x_2]. The partial dependence plot for variable x_1 is model m is function p_1 defined as x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)]. This can be approximated, using some dataset using \widehat{p}_1(x_1)=\frac{1}{n}\sum_{i=1}^n m(x_1,x_{2,i})My concern here what the interpretation of that plot when there are some (strongly) correlated covariates. Let us generate some dataset to start with

n=1000
library(mnormt)
r=.7
set.seed(1234)
X = rmnorm(n,mean = c(0,0),varcov = matrix(c(1,r,r,1),2,2))
Y = 1+X[,1]-2*X[,2]+rnorm(n)/2
df = data.frame(Y=Y,X1=X[,1],X2=X[,2])

As we can see, the true model is here is y_i=\beta_0+\beta_1 x_{1,i}+\beta_2x_{2,i}+\varepsilon_i where \beta_1 =1 but the two variables are positively correlated, and the second one has a strong negative impact. Note that here

reg = lm(Y~.,data=df)
summary(reg)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.01414    0.01601   63.35   <2e-16 ***
X1           1.02268    0.02305   44.37   <2e-16 ***
X2          -2.03248    0.02342  -86.80   <2e-16 ***

If we estimate a wrongly specified model y_i=b_0+b_1 x_{1,i}+\eta_i, we would get

reg1 = lm(Y~X1,data=df)
summary(reg1)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.03522    0.04680  22.121   <2e-16 ***
X1          -0.44148    0.04591  -9.616   <2e-16 ***

Thus, on the proper model, \widehat{\beta}_1\sim+1.02 while \widehat{b}_1\sim-0.44  on the mispecified model.

Now, let us look at the parial dependence plot of the good model, using standard R dedicated packages,

library(pdp) 
pdp::partial(reg, pred.var = "X1", plot = TRUE,
              plot.engine = "ggplot2")

which is the linear line y=1+x, that corresponds to y=\beta_0+\beta_1x.

library(DALEX)
plot(DALEX::single_variable(DALEX::explain(reg,
data=df),variable = "X1",type = "pdp"))

which corresponds to the previous graph. Here, it is also possible to creaste our own function to compute that partial dependence plot,

pdp1 = function(x1){
  nd = data.frame(X1=x1,X2=df$X2)
  mean(predict(reg,newdata=nd))
}

that will be the straight line below (the dotted line is the theoretical one y=1+x,

vx=seq(-3.5,3.5,length=101)
vpdp1 = Vectorize(pdp1)(vx)
plot(vx,vpdp1,type="l")
abline(a=1,b=1,lty=2)

which is very different from the univariate regression on x_1

abline(reg1,col="red")

Actually, the later is very consistent with a local regression, only on x_1

library(locfit)
lines(locfit(Y~X1,data=df),col="blue")

Now, to get back to the definition of the partial dependence plot, x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)], in the context of correlated variable, I was wondering if it would not make more sense to consider some local version actually, something like x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2|X_1}}[m(x_1,X_2)]. My intuition was that, somehow, it did not make any sense to consider any X_2 while X_1 was fixed (and equal to x_1). But it would make more sense actually to look at more valid X_2‘s given the value of X_1. And a natural estimate could be some k neareast-neighbors, i.e. \tilde{p}_1(x_1)=\frac{1}{k}\sum_{i\in\mathcal{V}_k(x)}^n m(x_1,x_{2,i}) where \mathcal{V}_k(x) is the set of indices of the k x_i‘s that are the closest to x, i.e.

lpdp1 = function(x1){
  nd = data.frame(X1=x1,X2=df$X2)
  idx = rank(abs(df$X1-x1))
  mean(predict(reg,newdata=nd[idx<50,]))
}
vlpdp1 = Vectorize(lpdp1)(vx)
lines(vx,vlpdp1,col="darkgreen",lwd=2)

Surprisingly (?), this local partial dependence plot gives a curve that corresponds to the simple regression…

IRE (HEC) – CREEI Seminar, HEC Montréal, on Monday

This Monday, I will give a talk at the IRE (HEC) – CREEI Seminar, at HEC Montréal. I will have one hour to present our work, with Olivier Cabrignac and Ewen Gallic, on “Modeling Joint Lives within Families

Les conséquences de la loi de Goodhart

Comme l’a énoncé Marilyn Strathern, la loi de Goodhart dit que « lorsqu’une mesure devient un objectif, elle cesse d’être une bonne mesure ». Les applications économiques sont nombreuses, mais cette loi permet aussi de comprendre les dangers des décisions algorithmiques, ou d’expliquer la difficulté à utiliser les données disponibles depuis le début de la pandémie de SARS-CoV-2 COVID-19.

La loi de Goodhart, l’évaluation des  politiques publiques et Facebook

« Lorsqu’une mesure devient un objectif, elle cesse d’être une bonne mesure » est la formulation la plus simple de la loi de Goodhart. Certains parlent aussi de loi de Campbell[i], Donald Campbell ayant affirmé « the more any quantitative social indicator is used for social decision-making, the more subject it will be to corruption pressures and the more apt it will be to distort and corrupt the social processes it is intended to monitor ». Cette corruption des métriques, ou des statistiques, a été observée dans de nombreux domaines, en particulier en santé, en matière de justice, ou en éducation.

Les budgets de certaines écoles étaient directement fonction des résultats obtenus à certains tests. On comprend que les enseignants sont fortement incités à commencer à préparer à l’examen au lieu d’enseigner de manière généraliste, voire à retirer discrètement certains enfants qui risqueraient d’échouer, comme le raconte Strauss (2015). En fixant le taux de réussite comme un objectif à maximiser, cette mesure ne signifie plus rien car elle a induit des comportements incitatifs qui ne reflètent plus la qualité de la formation. . Ce reproche peut également être fait à de nombreuses associations actuarielles à travers le monde, qui ont mis en place des ‘examens professionnels’, que peuvent passer les étudiants ayant eu une formation en actuariat. De nombreuses formations universitaires ont alors délaissé un enseignement généraliste pour offrir, à la place, des préparations à ces examens. Les examens ne sont plus une ‘mesure’ des connaissances des étudiants, ces derniers ayant de plus en plus de mal à sortir du cadre très scolaire des exercices très formatés.

Dans le domaine de la santé, aux États-Unis, Poku (2015) note qu’à partir de 2012, en vertu de la loi ‘Affordable Care Act’, Medicare a commencé à imposer des sanctions financières aux hôpitaux présentant des taux de réadmission de 30 jours « plus élevés que prévu ». Dès lors, le taux moyen de réadmission à l’hôpital dans les 30 jours pour les bénéficiaires de l’assurance fee-for-service a diminué. Est-ce dû à une amélioration des efforts de transition et de coordination des soins menées par les hôpitaux, ou bien faut-il relier cette baisse à une hausse des séjours « en observation » pendant la même période ? Car bien souvent, fixer un objectif sur la base d’une mesure précise (ici le taux de réadmission de 30 jours) rend cette grandeur totalement inutilisable pour quantifier le risque de retomber malade, mais a aussi une influence directe sur d’autres grandeurs (ici le nombre de séjours « en observation »), qu’il devient alors difficile de suivre sur la durée.

Sur internet, on demande de plus en plus à des algorithmes de trier du contenu, de juger le caractère diffamatoire ou raciste des tweets, de voir si une vidéo relève du deepfake, de mettre un score de fiabilité d’un compte Facebook, etc. Et nombreux sont ceux qui voudraient savoir comment ces scores sont créés. Malheureusement, comme le notait Dwoskin (2018) « not knowing how [Facebook is] judging us is what makes us uncomfortable. But the irony is that they can’t tell us how they are judging us — because if they do, the algorithms that they built will be gamed »[ii], exactement comme le suppose la loi de Goodhart.

La critique de Lucas et le thermostat de Friedman

Au début des années 1970, Robert Lucas expliquait que les décideurs économiques devaient éviter de se baser « naïvement » sur des statistiques passées pour prédire le comportement futur des agents, « given that the structure of an econometric model consists of optimal decision rules of economic agents, and that optimal decision rules vary systematically with changes in the structure of series relevant to the decision maker, it follows that any change in policy will systematically alter the structure of econometric models »[iii]. L’idée sous-jacente était simplement que les agents s’adaptent aux signaux qu’ils reçoivent.

Comme le dira Charles Goodhart quelques années plus tard « as the statistical relationships derived from the past depended on the particular kind of policy aim pursued by the authorities over the period considered, there would be no guarantee of their exact continuation in the future, should that policy be altered ». Ou formulé autrement, toute relation statistique observée (on pourra penser à une corrélation forte entre deux variables) aura tendance à disparaître une fois qu’une pression est exercée sur elle à des fins de contrôle. En fait, Charles Goodhart va plus loin que Robert Lucas dans Goodhart (1975), laissant entendre que dans de nombreux cas, les agents vont modifier leur comportement à leur avantage, même si cela se fait au détriment d’un éventuel bien-être collectif (on pourra repenser aux exemples de l’éducation, ou de la santé). Cette loi a été élaborée après avoir observé comment le gouvernement de Margaret Thatcher, dans les années 1980, a ciblé l’offre de monnaie pour contrôler l’inflation, mais a ensuite constaté que les agrégats monétaires avaient perdu leur relation précédemment forte avec l’inflation. L’inflation a échappé à tout contrôle, même lorsque le gouvernement exerçait une pression forte sur la masse monétaire.

Friedman (2003) avait utilisé l’analogie du thermostat pour expliquer le problème : la banque centrale est le thermostat de l’économie. Le décideur dispose d’une information ( ) lui permettant d’agir sur une variable de contrôle ( ), pour s’assurer qu’une variable ( ) soit proche de la valeur cible ( ). Et les erreurs de prévision ( ) se doivent d’être décorrélées de l’information  et du contrôle , si les anticipations sont rationnelles, ce qui peut sembler paradoxal. Pour reprendre l’image de Farrell (2012), imaginons un conducteur qui circule en voiture sur une route très vallonnée, à vitesse constante. Avec des hauts et des bas, cela signifie que le conducteur dose parfaitement ses accélérations et ses freins pour contrôler la vitesse. Pourtant, vu de loin, on serait tenté de dire que le niveau d’accélération n’a aucun impact sur la vitesse, et que si on faisait une régression de la vitesse sur l’accélération de la voiture, la corrélation entre les deux serait nulle, autrement dit, ici, accélérer et décélérer n’a aucune influence sur la vitesse…

Au-delà du contrôle

Le problème évoqué par la loi de Goodhart se retrouve aussi en modélisation statistique, et en apprentissage machine. Dans ce dernier cas, la fonction objectif que l’on optimise est en lien avec le problème que l’on cherche à résoudre : dans les algorithmes prédictifs de régression, on va comparer la prévision avec la réalisation, et pourra prendre la somme des carrés des erreurs ; pour un problème de classification ou de labellisation (fraude/pas fraude), on comptera le nombre d’erreur de classification (avec des fonctions de coût potentiellement différentes si les deux types d’erreurs ont des impacts différents). Mais bien souvent, on ne cherche pas un modèle parfait, sans erreur, on veut un modèle qui prédira bien sur des nouvelles données ! Aussi, on va éviter de juger des qualités prédictives d’un modèle sur les données qui ont été utilisées pour construire le modèle. On va alors utiliser une partie des données pour construire le modèle, et une autre pour juger de son pouvoir prédictif, et voir à partir de quel moment le modèle commence à modéliser du bruit, au lieu de chercher un lien fort entre les variables explicatives, et la variable d’intérêt. Cette approche, généralisée avec la notion de « validation croisée », permet de séparer l’objectif de la mesure.

En programmation dynamique, la loi de Goodhart est également bien connue (même si c’est parfois avec un autre nom). Dans les systèmes dynamiques, l’agent s’intéresse à une grandeur (x_t), par exemple un stock de biens qu’il vendra, et cherchera à maximiser une fonction de la forme f(x_1,x_2,\cdots,x_T), à partir d’une valeur initiale donnée, x_0, par exemple la somme (escomptée ou non) de tous les x_t, ou bien peut être uniquement la valeur terminale x_T. La dynamique de (x_t) dépend d’une variable de contrôle, (u_t), que l’agent peut choisir, en sachant que x_{t+1} dépendra directement de u_t, et éventuellement d’autres grandeurs, comme x_t. Bellman (1957) a posé les bases mathématiques de la résolution de ce genre de problème, que l’on retrouve généralisés que l’on retrouve dans l’apprentissage par renforcement[iv], où l’agent devra explorer, tenter différents contrôles, afin d’apprendre la manière dont u_t  va influencer x_{t+1}. Un exemple récent pourrait être le contrôle d’une pandémie, où x_t serait le nombre de personne infectées, ou le nombre de décès, et u_t est un levier de contrôle, comme le nombre de tests proposés, ou le nombre de personnes autorisés à aller travailler. Naturellement, si x_t mesure le nombre de personnes testées positives à la date t, une variable de contrôle permettant de faire baisser facilement x_t est le nombre de tests effectués, mais cela ne va en rien ralentir la propagation de l’épidémie (ce qui semble l’objectif naturel). En réalité, comme le notait l’immunologue Anthony Fauci, directeur depuis 1984 de l’Institut National des Allergies et Maladies Infectieuses aux États-Unis, « if it looks like you’re overreacting, your’re probably doing the right thing » (cité dans Budryk (2020)), compte tenu des effets de rétroaction.

De l’overdose de quantification et la maîtrise des métriques

Avant de revenir sur la pandémie de 2020, notons que les métriques sont souvent introduites dans un souci de transparence, de semblant d’exigence démocratique, comme la traduction ou l’expression concrète d’un objectif collectif, mais aussi souvent individuel. Mais toute métrique, toute statistique, cache souvent une réalité bien plus complexe. L’exemple des statistiques du chômage a probablement été un des plus étudiés, comme le montre Desrosières (2008). La série statistique conjoncturelle du chômage, publiée par l’INSEE, était très attendue par les responsables politiques et par la presse, au point de devenir l’objectif affiché de plusieurs gouvernements : « faire baisser les chiffres du chômage », comme le rappelle Errard (2015) par exemple. En voulant donner l’illusion de contrôler (et de faire baisser) le chômage, la pression a été mise pour que les conseillers de Pôle Emploi pour augmenter les radiations, pour proposer des formules pour les jeunes, pour encourager à prendre deux mi-temps plutôt qu’un temps plein. Une fois compris comment la métrique cible était calculée, il était alors facile de la faire baisser artificiellement. Comme l’écrivait Desrosières (2000), « les indicateurs quantitatifs rétroagissent sur les acteurs quantifiés », comme on l’évoquait dans la formalisation mathématique du contrôle optimal dynamique. Mais plus que l’économiste (ou l’économètre) Charles Goodhart, c’est surtout Donald Campbell, qui a cherché à comprendre la manière dont les métriques faussent le comportement, et amènent les participants à les exploiter. Et récemment, Bruno & Didier (2013), ou Muller (2018) montrent comment éviter que la loi de Goodhart ne s’applique.

Car comme le notait déjà Charles Goodhart, il arrive que l’optimisation explicite d’un système utilisant une métrique rende finalement la métrique inutilisable, car au final plus du tout corrélée avec l’objectif. C’est le cas de bon nombre de système de punitions et de récompenses, qui visent à créer des incitations. On peut penser à l’assiduité en classe, par exemple, ou des absences sont punies afin d’inciter les élèves à travailler. En mettant en place de telles mesures, la fréquentation en classe augmente, mais les élèves ne travaillent pas pour autant. Et ceci même si une relation causale existe, entre la mesure et l’objectif, et pas une simple corrélation. Comme dans l’exemple de Facebook, garder une métrique secrète est une solution facile, une autre étant d’utiliser plusieurs métriques.

La loi de Goodhart et la crise sanitaire actuelle

Un des objectifs, régulièrement martelé depuis mars 2020, est qu’il ne faut pas saturer les systèmes de santés, dans tous les pays, le fameux « flatten the curve », évoqué par Ferguson et al. (2020). Il semblait essentiel de garantir, à tout prix, que les hôpitaux ne soient pas débordés. Au printemps 2020, les chaînes télévisées d’information donnaient, en temps continu, le nombre de personnes en soin intensifs, et le nombre de morts dans les hôpitaux, mesures que l’on retrouvera ensuite sous forme de graphiques, mis à jour toutes les semaines, voire tous les soirs, sur des sites internet dédiés. En cette période de crise, au plus fort de la saturation des hôpitaux, le N.H.S. en Angleterre avait demandé à chaque hôpital d’estimer ses capacités de lits, afin de réallouer les ressources globalement. Annoncer que peu de lits étaient disponibles était la meilleure stratégie pour obtenir davantage de financement. On peut alors s’interroger sur le niveau de saturation réelle du système, chaque hôpital ayant compris la règle, et manipulant la mesure à sa guise. Et tout aussi préoccupant, alors que les gouvernements se concentraient sur les hôpitaux (fournissant les données officielles utilisées pour construire la plupart des indicateurs), les maisons de retraites ont connu des hécatombes désastreuses, qui ont mis beaucoup de temps à être quantifiées. Giles (2020) raconte qu’en Angleterre, certains médecins auraient demandé à leurs patients âgés de réfléchir sérieusement à la question de savoir s’ils voulaient vraiment aller à l’hôpital, et utiliser les services d’urgence, au risque de passer plusieurs semaines isolés, coupés de leur famille.

Les statistiques sur le nombre de personnes (officiellement) positives n’ont cessé de dérouter les statisticiens, car aisément manipulables. On a tous en mémoire les déclarations de Donald Trump au début de l’été 2020, repris par exemple par Sheth (2020), qui affirmait que pour faire baisser le nombre de personnes ‘positives’ il suffisait de tester moins. Au début de la pandémie, un objectif clairement énoncé était de détecter les positifs asymptomatiques, et donc un ciblage des tests s’imposait. Le taux de positif était alors le signe que le ciblage des tests fonctionnait. En revanche, pour suivre l’évolution de la pandémie, il était indispensable d’effectuer les tests de manière aussi aléatoire que possible.

La crise engendrée par la pandémie SARS-CoV-2 COVID-19, avec sa quantification à outrance, avec ses statistiques mises à jour en temps réel, nous a rappelé les dangers de la loi de Goodhart. Comme le notait Laroussie (2021), le flot de statistiques a aussi permis à bon nombre d’entre nous de tenter l’exercice, d’essayer de prévoir l’évolution future des courbes, mais aussi de s’interroger sur la fiabilité des données, et de leur construction. Suivre le nombre de personnes supposées positives sans comprendre qui était testé, avec quel type de test, n’avait alors aucun sens. La dynamique des courbes était alors elle-même impactée par une boucle de rétroaction, découlant de décisions de décideurs politiques, qui avaient décidé de tester moins de personnes âgées lorsqu’il a fallu reprendre le chemin du travail, par exemple. Comment prendre des décisions avisées de matière de politique publique dans ces conditions ? C’est finalement la question profonde que pose la loi de Goodhart, nous rappelant aussi que les décideurs doivent apprendre à faire la différence entre l’esprit de la loi et la lettre de la loi – les routes de l’enfer étant pavées de bonnes intention – en gardant un esprit de mesure.

Références

Bellman, R. (1957). Dynamic Programming. Princeton University Press.

Bruno, I. & Didier, E. (2013). Benchmarking. L’État sous pression statistique. Paris, La Découverte.

Budryk, Z. (2020). ‘If it looks like you’re overreacting, you’re probably doing the right thing‘ The Hill, 15 mars 2020.

Campbell, D. T. (1975) Assessing the impact of planned social change. In G. M. Lyons

(ed.), Social Research and Public Policies: The Dartmouth/OECD Conference (pp. 3–45). Hanover, NH: Public Affairs Center

Charpentier, A., Elie, R. & Remlinger, C. (2020). Reinforcement Learning in Economics and Finance. arXiv:2003.10014

Daston, L. (2010). Why statistics tend not only to describe the world but to change it. The London Review of Books, 22:8.

Desrosières, A. (2000). La Politique des grands nombres : Histoire de la raison statistique. La Découverte.

Desrosières, A. (2008). Gouverner par les nombres. Presses de l’École des Mines.

Dwoskin, E. (2018). Facebook is rating the trustworthiness of its users on a scale from zero to one. Washington Post, 21 aout 2018,

Errard, G. (2015). Le contrôle des chômeurs peut-il faire baisser le chômage ?. Le Figaro, 26 août 2015,

Farrell, H. (2012). Milton Friedman’s Thermostat. Monkey Cage, 31 juillet 2012.

Ferguson, N. et al. (2020). Impact of non-pharmaceutical interventions to reduce covid-19 mortality and healthcare demand. Imperial College COVID-19 Response Team 9.

Friedman, M. (2003). The Fed’s Thermostat. The Wall Street Journal, 19 août 2003.

Giles, C. (2020). Goodhart’s law comes back to haunt the UK’s Covid strategy. Financial Times, 14 mai 2020,

Goodhart, C.A.E. (1975) Problems of monetary management: The UK experience. Papers in Monetary Economics, Volume I. Sydney: Reserve Bank of Australia.

Laroussie, D. (2021). Covid-19 : ces modélisateurs qui anticipent la pandémie. Le Monde, 5 janvier 2021,

Muller, J. Z. (2018). The tyranny of metrics. Princeton University Press.

Rodamar, J. (2018) There ought to be a law! Campbell versus Goodhart. Significance, 15:6.

Sheth, S. (2020). Trump says that ‘if we stop testing right now, we’d have very few cases’ of the coronavirus. Business Insider, 15 juin 2020

Strauss, V. (2015). How and why convicted Atlanta teachers cheated on standardized tests. The Washington Post, 1 avril 2015,

[i] Rodamar (2018) revient sur la comparaison entre les deux publications, Goodhart (1975) et Campbell (1975), qui énoncent le même principe, dans des contextes très différents.

[ii] « ne pas savoir comment [Facebook] nous juge, c’est ce qui nous met mal à l’aise. Mais l’ironie est qu’ils ne peuvent pas nous dire comment ils nous jugent – parce que s’ils le font, les algorithmes qu’ils ont construits seront alors manipulés » (traduction de l’auteur)

[iii] « étant donné que la structure d’un modèle économétrique consiste en des règles de décision optimales des agents économiques, et que les règles de décision optimales varient systématiquement en fonction des changements dans la structure des séries pertinentes pour le décideur, il s’ensuit que tout changement de politique modifiera systématiquement la structure des modèles économétriques » (traduction de l’auteur).

[iv] Ou reinforcement learning, comme décrit dans Charpentier et al.  (2020).

Lilliefors, Kolmogorov-Smirnov and cross-validation

In statistics, Kolmogorov–Smirnov test is a popular procedure to test, from a sample \{x_1,\cdots,x_n\} is drawn from a distribution F, or usually F_{\theta_0}, where F_{\theta} is some parametric distribution. For instance, we can test H_0:X_i\sim\mathcal{N(0,1)} (where \theta_0=(\mu_0,\sigma_0^2)=(0,1)) using that test. More specifically, I wanted to discuss today p-values. Given n let us draw \mathcal{N}(0,1) samples of size n, and compute the p-values of Kolmogorov–Smirnov tests

n=300
p = rep(NA,1e5)
for(s in 1:1e5){
X = rnorm(n,0,1)
p[s] = ks.test(X,"pnorm",0,1)$p.value
}

We can visualise the distribution of the p-values below (I added some Beta distribution fit here)

library(fitdistrplus)
fit.dist = fitdist(p,"beta")
hist(p,probability = TRUE,main="",xlab="",ylab="")
vu = seq(0,1,by=.01)
vv = dbeta(vu,shape1 = fit.dist$estimate[1], shape2 = fit.dist$estimate[2])
lines(vu,vv,col="dark red", lwd=2)

It looks like it is quite uniform (theoretically, the p-value is uniform). More specifically, the p-value was lower than 5% in 5% of the samples

[note: here I compute ‘mean(p<=.05)’ but I have some trouble with the ‘<‘ and ‘>’ symbols, as always]

mean(p&lt;=.05)
[1] 0.0479

i.e. we wrongly reject H_0:X_i\sim\mathcal{N(0,1)} is 5% of the samples.

As discussed previously on the blog, in many cases, we do care about the distribution, and not really the parameters, so we wish to test something like H_0:X_i\sim\mathcal{N(\mu,\sigma^2)}, for some \mu and \sigma^2. Therefore, a natural idea can be to test H_0:X_i\sim\mathcal{N(\hat\mu,\hat\sigma^2)}, for some estimates of \mu and \sigma^2. That’s the idea of Lilliefors test. More specifically, Lilliefors test suggests to use , Kolmogorov–Smirnov statistics, but corrects the p-value. Indeed, if we draw many samples, and use Kolmogorov–Smirnov statistics and its classical p-value to test for H_0:X_i\sim\mathcal{N(\hat\mu,\hat\sigma^2)},

n=300
p = rep(NA,1e5)
for(s in 1:1e5){
X = rnorm(n,0,1)
p[s] = ks.test(X,"pnorm",mean(X),sd(X))$p.value
}

we see clearly that the distribution of p-values is no longer uniform

fit.dist = fitdist(p,"beta")
hist(p,probability = TRUE,main="",xlab="",ylab="")
vu = seq(0,1,by=.01)
vv = dbeta(vu,shape1 = fit.dist$estimate[1], shape2 = fit.dist$estimate[2])
lines(vu,vv,col="dark red", lwd=2)

More specifically, if x_i‘s are actually drawn from some Gaussian distribution, there are no chance to reject H_0, the p-value being almost never below 5%

mean(p&lt;=.05)
[1] 0.00012

Usually, to interpret that result, the heuristics is that \hat\mu and \hat\sigma^2 are both based on the sample, while previously 0 and 1 where based on some prior knowledge. Somehow, it reminded me on the classical problem when mention when we introduce cross-validation, which is Goodhart’s law

When a measure becomes a target, it ceases to be a good measure

i.e. we cannot assess goodness of fit using the same data as the ones used to estimate parameters. So here, why not use some hold-out (or cross-validation) procedure : split the dataset in two parts, \{x_1,\cdots,x_k\} (with k<n) to estimate parameters \mu and \sigma^2 and then use \{x_{k+1},\cdots,x_n\} and Kolmogorov–Smirnov statistics on it to test if x_i‘s are drawn from some Gaussian distribution. More precisely, will the p-value computed using the standard Kolmogorov–Smirnov procedure be ok here. Here, I tried two scenarios, k/n being either 1/3 or 2/3,

p = matrix(NA,1e5,4)
for(s in 1:1e5){
X = rnorm(n,0,1)
p[s,1] = ks.test(X,"pnorm",0,1)$p.value
p[s,2] = ks.test(X,"pnorm",mean(X),sd(X))$p.value
p[s,3] = ks.test(X[1:200],"pnorm",mean(X[201:300]),sd(X[201:300]))$p.value
p[s,4] = ks.test(X[201:300],"pnorm",mean(X[1:200]),sd(X[1:200]))$p.value
}

Again, we can visualize the distributions of p-values,  in the case where 1/3 of the data is used to estimate \mu and \sigma^2, and 2/3 of the data is used to test

fit.dist = fitdist(p[,3],"beta")
hist(p[,3],probability = TRUE,main="",xlab="",ylab="")
vu=seq(0,1,by=.01)
vv=dbeta(vu,shape1 = fit.dist$estimate[1], shape2 = fit.dist$estimate[2])
lines(vu,vv,col="dark red", lwd=2)


and in the case where 2/3 of the data is used to estimate \mu and \sigma^2, and 1/3 of the data is used to test

fit.dist = fitdist(p[,4],"beta")
hist(p[,4],probability = TRUE,main="",xlab="",ylab="")
vu=seq(0,1,by=.01)
vv=dbeta(vu,shape1 = fit.dist$estimate[1], shape2 = fit.dist$estimate[2])
lines(vu,vv,col="dark red", lwd=2)


Observe here that we (wrongly) reject too frequently H_0, since the p-values are  below 5% in 25% of the scenarios, in the first case (less data used to estimate), and 9% of the scenarios, in the second case (less data used to test)

mean(p[,3]&lt;=.05)
[1] 0.24168
mean(p[,4]&lt;=.05)
[1] 0.09334

We can actually compute that probability as a function of k/n

n=300
p = matrix(NA,1e4,99)
for(s in 1:1e4){
  X = rnorm(n,0,1)
  KS = function(p) ks.test(X[1:(p*n)],"pnorm",mean(X[(p*n+1):n]),sd(X[(p*n+1):n]))$p.value
  p[s,] = Vectorize(KS)((1:99)/100)
}

The evolution of the probability is the following

prob5pc = apply(p,2,function(x) mean(x&lt;=.05))
plot((1:99)/100,prob5pc)

so, it looks like we can use some sort of hold-out procedure to test for H_0:X_i\sim\mathcal{N(\mu,\sigma^2)}, for some \mu and \sigma^2, using Kolmogorov–Smirnov test with \mu=\hat\mu and \sigma^2=\hat\sigma^2 but the proportion of data used to estimate those quantities should be (much) larger that the one used to compute the statistics. Otherwise, we clearly reject too frequently H_0.

Webinaire sur l’assurance collaborative et la théorie des grahs

Le mardi 28 janvier, avec Christian Yann Robert, on donnera un exposé à Paris sur le thème “Assurance collaborative: Théorie des Graphes, Machine Learning et Actuariat“, à l’Institut Louis Bachelier et l’Institut des Actuaires, à Paris

L’assurance collaborative est une nouvelle forme d’assurance et de relation entre assureur et assurés qui met à profit les nouvelles technologies pour une plus grande implication des assurés dans le processus et le modèle économique. Par la transparence elle vise à instaurer une relation de confiance forte entre les acteurs. Elle peut prendre des formes différentes, allant de la simple redistribution des profits jusqu’à la gestion en peer to peer par les assurés, rassemblés en communautés d’intérêt. Plusieurs tentatives ont été faites depuis 5 ou 6 ans dans différents pays: certaines en Europe ont échoué, d’autres se rôdent.
Aux Etats-unis, Lemonade est une vraie réussite. La répartition de la prime est de 20% pour la gestion et la rémunération de l’entreprise, 40% pour les sinistres (avec une redistribution en cas de bénéfice technique) et 40% pour la réassurance des gros sinistres. Cette nouvelle forme de mutualisation fait émerger de nouvelles problématiques tarifaires et de provisionnement, selon le mode de constitution des communautés d’assurés. La théorie des graphes aide à comprendre la capacité de mutualisation selon la structure des communautés.