Lissage

Courbe lissée d'un ensemble de points

On suppose que l'on dispose d'un ensemble de points de mesures bruitées. Pas de modèle statistique simple des données, mais on veut juste en tirer une courbe.

N <- 1000 #On va travailler sur 1000 points 
x <- runif(N,min=0,max = 10) # valeurs de x dans [0;1]
y <- x^2 + sqrt(x) + 10*cos(x) + 5*rnorm(N) # La formule que l'on ne connait pas
plot(x,y,col=rgb(0,100,0,70,maxColorValue=255),pch=16) 

Ce qui donne un truc du genre de : (les points sont affichés avec une opacité de 70%, c'est bien pratique pour visualiser des ensembles de points un peu compacts).

Ensemble de points

L'astuce est de faire une convolution par une gaussienne. C'est facile s'adapte en dimension plus grande (au prix de pas mal de calculs en plus). Evidemment il faudra choisir la largeur de bande.

#Le noyau de la convolution
K<-function(x,s){
1/(sqrt(2*pi)*s)*exp(-x^2/s)
}

         
sigma <- 1
support = seq(min(x),max(x),by = 0.01 )

#Lissage. L'on sait que les bords vont poser des problèmes...
conv <- c()
for (i in 1:length(support)){
	buf <- K( x - support[i], sigma )
	conv[i] <- sum(y * buf) / sum(buf)
	}

#Voila c'est déjà fini.
#Affichage de la courbe lissée	
plot(x,y,col=rgb(0,100,0,70,maxColorValue=255),pch=16) 	
lines(support,conv,col='red',lwd=2)	

#idéalement
X <- sort(x)
Y <- X^2 + sqrt(X) + 10*cos(X) 
lines(X,Y,col='blue',lwd=2)
legend("topleft",legend=c('Courbe idéale','Lissage'), col = c('blue','red'), lwd =  2)

Et un résultat possible :

lissage

Avec un sigma de 0.15 c'est encore mieux :

lissage

Pour le fun, on recommence en 2D avec cette surface comme objectif (prise dans la doc de R)

x <- seq(-10, 10, length= 60)
y <- x
f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
z <- outer(x, y, f)
z[is.na(z)] <- 1
op <- par(bg = "white")
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Sinc( r )"
) -> res
round(res,3)

objectif

On bruite les points

Z <- z + matrix(rnorm(60*60,sd=1),60,60)
persp(x, y, Z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Sinc( r ) + bruit"
) -> res
round(res,3)

points bruités

Et on lisse (en utilisant le même support qu'au départ). Attention quand même içi on utilise un support trop grand pour la gaussienne. En traitement d'image on appelle ce genre de traitement un filtrage. De plus dans ce cas particulier la convolution 2D peut être ramenée à 2 convolutions 1D.

K<-function(x,y,s=1){
1/(sqrt(2*pi)*s)*exp(-(x^2+y^2)/s)
}

conv <- 0 * Z;
for (i in 1:length(x)){
	for (j in 1:length(y)){
		cum <- 0 #A cause des bords
		for (ii in 1:length(x)){
			for (jj in 1:length(y)){
				w <- K( x[ii] - x[i], y[jj]- y[j] )
				cum <- cum + 	w
				conv[i,j] <- conv[i,j] + Z[ii,jj]*w
			}
		}
		conv[i,j]<-conv[i,j]/cum
	}
}

#affichage
persp(x, y, conv, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Lissage"
) -> res
round(res,3)

Lissage

Remarque : on peut faire beaucoup mieux en utilisant le fait que la fonction possède une invariance de rotation (et est donc en fait une fonction d'une seule variable). C'est évidemment une hyppothèse supplémentaire fausse dans le cas général. Je le laisse en exercice.

Design selector