La méthode de Whittaker-Henderson est une méthode de lissage non paramétrique ; combinant un critère de fidélité aux données et un critère de régularité, elle permet de lisser simplement des surfaces et est à ce titre utile notamment en arrêt de travail.

Références

Support de cours sur les méthodes de lissage et d'ajustements

Mise en oeuvre

La fonction ci-dessous permet de réaliser un lissage en dimension 2 :

W_H<-function(Tbrut,Poids,Ordre_vertical,Ordre_horizontal,alpha,beta){

#Déclaration de la taille des variables
q<<-ncol(Tbrut)
p<<-nrow(Tbrut)

#Définition des matrices
U=matrix(0,q*p)
v=matrix(0,Ordre_vertical+1)
h=matrix(0,Ordre_horizontal+1)
x=matrix(1:p,p,1)
y=matrix(1:q,q,1)
Kv=matrix(0,(p-Ordre_vertical)*q,q*p)
Kh=matrix(0,p*(q-Ordre_horizontal),q*p)
M=matrix(0,q*p,q*p)
W=matrix(0,q*p,q*p)
qlisse=matrix(0,p,q)
Tlisse=matrix(0,p,q)
Ecart=matrix(0,p,q)

#Transformation de la matrice de taux bruts en vecteur et construction de la matrice des poids
for (j in 1:q) {
for (i in 1:p) {
U[(i-1)*q+j,1]=Tbrut[i,j]
W[q*(i-1)+j,q*(i-1)+j]=Poids[i,j]
}
}

#Construction matrice kv
for (k in 0:Ordre_vertical) {
v[(k+1),1]=(-1)^(Ordre_vertical-k)*factorial(Ordre_vertical)/(factorial(k)*factorial(Ordre_vertical-k))
}

for (j in 1:q) {
for (z in 1:(p-Ordre_vertical)) {
for (i in 1:(Ordre_vertical+1)) {
Kv[z+(j-1)*(p-Ordre_vertical),j+(q)*(i-1)+(z-1)*(q)]=v[i,1]
}
}
}

#Construction matrice Kh
for (k in 0:Ordre_horizontal) {
h[(k+1),1]=(-1)^(Ordre_horizontal-k)*factorial(Ordre_horizontal)/(factorial(Ordre_horizontal-k)*factorial(k))
}
for (i in 1:p) {
for (j in 1:(q-Ordre_horizontal)) {
for (z in 1:(Ordre_horizontal+1)) {
Kh[j+(i-1)*(q-Ordre_horizontal),z+(j-1)+(i-1)*(q)]=h[z,1]
}
}
}

#Calcul des taux lissés
M=W+alpha*t(Kv)%*%Kv+beta*t(Kh)%*%Kh
qlisse=solve(M)%*%W%*%U
for (j in 1:q) {
for (i in 1:p) {
Tlisse[i,j]=qlisse[(i-1)*q+j,1]
Ecart[i,j]=Tlisse[i,j]/Tbrut[i,j]
}
}

#Stockage des résultats
Tlisse<<-Tlisse
Ecart<<-Ecart

#Représentation graphique
par(mfrow=c(1,3))
persp(x,y,Tbrut,theta=45,col=brewer.pal(9,"RdYlGn"),xlab="Age",ylab="Ancienneté",main="Taux bruts",adj=0.5,font=2)
persp(x,y,Tlisse,theta=45,col=brewer.pal(9,"RdYlGn"),xlab="Age",ylab="Ancienneté",main="Taux lissés",adj=0.5,font=2)
persp(x,y,Ecart,theta=90,phi=22.5,expand=1,col=brewer.pal(9,"RdYlGn"),xlab="Age",ylab="Ancienneté",main="Taux lissés/Taux bruts",adj=0.5,font=2)

}

Ce code est fourni dans le fichier .r ci-joint.