Régression linéaire, méthode rapide

**progfou** · 12/02/2011, 10h07

Bonjour,

J'ai un ensemble de points (30 à 50) et je sais que je peux faire une régression linéaire pour déterminer une fonction affine.
Cependant, j'ai des contraintes de rapidité et de faible puissance (système embarqué).
Je cherche donc une méthode rapide pour faire cette estimation. Les moindres carrés sont couteux (en tout cas, je n'ai pas trouvé d'algorithme rapide).
J'ai eu une idée qui me paraissait assez simple et je voudrais avoir votre avis.
J'ai mes données dans un tableau, disons tab

avec N valeurs et l'indice commencant à 0.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
Trier tab
a1 = (tab[0] - tab[N])/(N-1)
a2 = (tab[1] - tab[N-1])/(N-3)
a = (a1+a2)/2

Evidemment, à faire sur plusieurs valeurs afin d'éliminer une partie du bruit (Avec mes données, il semblerait qu'avec 4 valeurs ça suffit).
Est-ce que ça estime effectivement mon coefficient a de y=a*x+b ?

Merci pour vos avis !

**Aleph69** · 12/02/2011, 12h28

Bonjour,

Est-ce que ça estime effectivement mon coefficient a de y=a*x+b ?

De quelle droite y parles-tu?

J'ai un ensemble de points (30 à 50) et je sais que je peux faire une régression linéaire pour déterminer une fonction affine.
Cependant, j'ai des contraintes de rapidité et de faible puissance (système embarqué).
Je cherche donc une méthode rapide pour faire cette estimation. Les moindres carrés sont couteux (en tout cas, je n'ai pas trouvé d'algorithme rapide).

Pour 30 à 50 points, le coût des moindres carrés est négligeable quand même. As-tu des contraintes de temps ou de complexité?

Pour information, les coefficients a et b de la droite y donnée par les moindres carrés sont définis explicitement ici :
http://fr.wikipedia.org/wiki/R%C3%A9..._lin%C3%A9aire

**b_reda31** · 13/02/2011, 10h24

Envoyé par progfou

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
Trier tab
a1 = (tab[0] - tab[N])/(N-1) 
a2 = (tab[1] - tab[N-1])/(N-3)
a = (a1+a2)/2

Si je comprend bien vous approximez la pente de la droite que vous recherchez par la moyenne des deux pente a1 et a2 :
a1 étant la pente de la droite passant par tab[0] et tab[N]. a2 la pente de la droite passant par tab[1]et tab[N-1]. Il me semble alors qu'il y a une erreur dans vos formules. Ne serait ce pas :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
Trier tab
a1 = (tab[0] - tab[N])/N 
a2 = (tab[1] - tab[N-1])/N-2
a = (a1+a2)/2

Ceci dit, ici vous réduisez l'information que porte l'ensemble de points (30 à 50 points) à celle de 4 points! Qu'en est il des autres points, ils ont aussi leur mot à dire non?
En supposant que chaque point dispose d'une certaine erreur d'estimation, si l'erreur dont dispose les quatre points tab[0],tab[1],tab[N-1],tab[N] est relativement importante, cela faussera complétement vos calculs il me semble.

**Aleph69** · 13/02/2011, 15h29

Bonjour,

une petite remarque au passage : la régression linéaire étant très sensible aux outliers, il vaudrait mieux voir du côté de la régression PLS.

**progfou** · 13/02/2011, 15h47

Merci pour vos réponses !
Effectivement l'idée était bien de récupérer seulement les extrèmes (par tri du tableau) afin de limiter les calculs.
Il faudrait que l'estimation soit faite en quelque dizaines de µs (une fois le tableau de points rempli), sur une machine d'une puissance relativement limitée (pour aller vite, additions, soustractions OK, multiplications et divisions à éviter).

Je regarde la PLS

**Aleph69** · 13/02/2011, 16h04

Bonjour,

plus tu vas réduire le nombre de points, plus tu vas détériorer la qualité de la régression. Pourquoi ne pas construire un modèle expert? Par exemple, une idée naïve : tu prends plusieurs nuages de points concrets et tu fais une belle régression en récupérant les coefficients a et b des droites y=ax+b pour chaque nuage. Ensuite, tu construits un modèle d'apprentissage pour prédire les valeurs de a et b. Le coût de l'apprentissage est très cher vis-àvis de tes contraintes mais est fait en amont, une fois pour toutes. En revanche, la prédiction statistique de a et b sera peu coûteuse pour les nouveaux nuages de points.

**souviron34** · 14/02/2011, 02h18

@progfou :

dans ta méthode, d'une part tu ne trouves pas vraiment (du tout) un coeff linéaire, d'autre part tu as non seulement un tri mais des divisions qui sont plus coûteuses que des multiplications....

Je sais que pseudocode avait mis en Java dans Contribuez une méthode de moindres carrés, et j'en avais mise une ici (En Fortran mais facilement adaptable) contenant une régrssion linéaire simple et efficace..

Sinon, en prenant ton idée de base, une manière plus précise et peu lourde serait, sans faire de tri initial :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
 
moy1 = 0.0
moy2 = 0.0
 
for ( i = 0 ; i < N/2 ; i++ )
      moy1 = moy1 + val[i] ;
 
for ( i = N/2 ; i < N ; i++ )
      moy2 = moy2 + val[i] ;
 
a = 2.0 * (moy2 - moy1) / N^2

Déjà, tu prends la moyenne sur la moitié de l'échantillon, donc c'est pas parfait, mais c'est très nettement mieux que sur 2 points...

Mais regarde quand même les moindres carrés...

**progfou** · 14/02/2011, 17h16

Envoyé par souviron34

@progfou :

dans ta méthode, d'une part tu ne trouves pas vraiment (du tout) un coeff linéaire, d'autre part tu as non seulement un tri mais des divisions qui sont plus coûteuses que des multiplications....

Je sais que pseudocode avait mis en Java dans Contribuez une méthode de moindres carrés, et j'en avais mise une ici (En Fortran mais facilement adaptable) contenant une régrssion linéaire simple et efficace..

Sinon, en prenant ton idée de base, une manière plus précise et peu lourde serait, sans faire de tri initial :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
 
moy1 = 0.0
moy2 = 0.0
 
for ( i = 0 ; i < N/2 ; i++ )
      moy1 = moy1 + val[i] ;
 
for ( i = N/2 ; i < N ; i++ )
      moy2 = moy2 + val[i] ;
 
a = 2.0 * (moy2 - moy1) / N^2

Déjà, tu prends la moyenne sur la moitié de l'échantillon, donc c'est pas parfait, mais c'est très nettement mieux que sur 2 points...

Mmh, intéressant, mais ça n'a pas l'air de fonctionner...
moy1 et moy2 sont en faite les sommes partielles (somme de la première moitié, et somme de la seconde) et non des moyennes.

Si je prends ces données là (mesurées) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Que j'applique l'algo, j'arrive à :
moy1 = 469774
moy2 = 474077

D'où : 2* (moy2 - moy1) / N^2 = 2*4303/(30^2) = 8606/900 = 9.56

Une courbe de régression donnerait un coefficient à 19.04.
On dirait qu'il y a un ratio de 2 entre les deux... normal ?

**souviron34** · 14/02/2011, 17h53

Envoyé par progfou

D'où : 2* (moy2 - moy1) / N^2 = 2*4303/(30^2) = 8606/900 = 9.56

Une courbe de régression donnerait un coefficient à 19.04.
On dirait qu'il y a un ratio de 2 entre les deux... normal ?

Sans doute....

Disons que j'avais divisé par N pour "égaliser"..

Mais si on considère que la moyenne moy1 s'applique sur l'intervalle 0-N/2, son abscisse devrait être à N/4, et à 3N/4 pour moy2..

Ma formule doit donc être fausse, et demander de diviser par N/2 et non par N..

D'où un facteur 2...

Et le point de départ est donc à b = moy1 - (moy2 - moy1)/2

**souviron34** · 14/02/2011, 18h16

Envoyé par progfou

moy1 et moy2 sont en faite les sommes partielles (somme de la première moitié, et somme de la seconde) et non des moyennes.

non, c'était des moyennes (voir ci-dessus le facteur 2/N)....

Pourquoi faire un calcul inutile ?

Au lieu de faire :

moy1 = moy1 / N/2
moy2 = moy2 / N/2

a = (moy2 - moy1) / N

(en reprenant la même erreur que j'avais faite),

on peut tout faire en une seule opération...

Mais comme on a besoin de moy1 (réel) pour le point de départ, on économise peu...

Donc si on a besoin de b vaut mieux calculer les moyennes réelles...

Mais si on a juste besoin de a, alors la formule

a = 4 * (moy2 - moy1) / N^2

avec moy1 et moy2 sommes partielles suffit...

**progfou** · 14/02/2011, 22h04

OK, ça me parait intéressant

.
Adapté avec un calcul virgule fixe, ça se passe bien

.
Il me reste un truc qui me paraît un peu compliqué : j'ai quelques points vraiment éloignés des autres (des outliers quoi). Ya moyen de m'en débarrasser facilement ?

**pseudocode** · 15/02/2011, 00h10

Envoyé par progfou

Les moindres carrés sont couteux (en tout cas, je n'ai pas trouvé d'algorithme rapide).

couteux, couteux...

a = (N.Sxy - Sx.Sy) / (N.Sxx - (Sx)²)

Sx = somme (x)
Sy = somme (y)
Sxy = somme (x.y)
Sxx = somme (x²)
N = nombre de points

Si ca prend trop de temps, il suffit de réduire le nombre de points.

**souviron34** · 15/02/2011, 00h23

Envoyé par progfou

Ya moyen de m'en débarrasser facilement ?

Facilement pas vraiment avec tes contraintes...

Soit un sigma soit un filtre...

Le plus simple à mettre en oeuvre est un sigma...

Ou effectuer d'abord une moyenne glissante, par exemple sur 3 (ou 5) points...

Mais ça fera des opérations...

Sinon une sorte de sigma non statistique....

Admettons que l'ordonnée de départ théorique soit ys, de fin ye... (calculé comme ci-dessus avec moy1 et moy2.)

Là il y a plusieurs choix :

un choix "hard-codé" dépendant de tes données : si tu sais qu'en général l'écart moyen est de 10% par exemple, tu précalcules 10% de la pente, et tu regardes si ton point correspond à l'ordonnée théorique +/- ce chiffre...

le choix plus statistique du vrai sigma... qui te demandera d'abord 3 boucles de plus (calculer la moyenne totale et la diff , puis repasser dans les points pour éliminer)..

Ou bien le faire par 1/2 comme pour les pentes...

Une moyenne glissante serait plus rapide à mettre en oeuvre et "aplatirait" la courbe, à condition que tu n'aies pas 3 outliers du même côté qui se suivent...

Mais nécessiterait pour faire bien un second tableau....

Donc je ne crois pas de solutions miracles de ce côté-là...

**progfou** · 15/02/2011, 08h36

Envoyé par pseudocode

couteux, couteux...
[...]
Sxy = somme (x.y)
Sxx = somme (x²)

C'est surtout ces deux sommes là qui posent problème, tant du point de vue calculs que du point de vue stockage (S(x.y) peut être grand !).
Ceci dit, il est peut être possible de le faire de façon maline... ?
Typiquement, les x sont connus, ils vont de 1 à N. Donc S(x²) peut être pré-calculé.
Ensuite, S(x.y) pour n=1 à N c'est :
y1 + 2.y2 + 3*y3 + ... + N*yn. Il est pas possible d'avoir une décomposition qui limite le nombre de multiplications ?
Déjà, tous les n puissances de 2 sont des décalages.
3 c'est un décalage et une addition.
Bref, un moyen algorithmique (une astuce mathématique ?) pour essayer de se rapprocher d'un nombre de multiplications le plus faible possible ?

**Aleph69** · 15/02/2011, 11h05

Bonjour,

si tu peux calculer les carrés des normes plus rapidement que les produits scalaires, alors tu peux exploiter les formules

||u+v||² = ||u||² + 2(u,v) + ||v||²,
||u-v||² = ||u||² - 2(u,v) + ||v||²,

pour deux vecteurs u et v de même longueur.

Sinon, pour la somme

y1 + 2.y2 + 3*y3 + ... + N*yn

ça se réécrit sous forme d'addition un peu à la façon Hörner :

(y1+y2+...+yn)+(y2+y3+...+yn)+(y3+y4+...+yn)+...+(y_{n-1}+yn)+yn,

soit de maniere equivalente

S1+S2+S3+...+S_{n-1}+Sn,

avec Si = somme des yj pour j allant de i à n.

EDIT : tu peux aussi exploiter les différences entre les Si. Tu construis une suite récurrente en posant y0=0 et Ti = Si-y_{i-1} pour i variant. Ca te permet de calculer la somme des yi une fois pour toute (i.e. S1) puis d'obtenir les TI en soustrayant un y à chaque fois. La formule devient
T1+T2+...+Tn.
Par contre, numériquement, tu dois perdre en stabilité mais pour N relativement petit ça devrait aller.

**pseudocode** · 15/02/2011, 11h53

Envoyé par progfou

Ensuite, S(x.y) pour n=1 à N c'est :
y1 + 2.y2 + 3*y3 + ... + N*yn. Il est pas possible d'avoir une décomposition qui limite le nombre de multiplications ?

En utilisant une sommation triangulaire ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
   y(n)
+  y(n) + y(n-1)
+  y(n) + y(n-1) + y(n-2)
+  y(n) + y(n-1) + y(n-2) + y(n-3)
+  ...
+  y(n) + y(n-1) + y(n-2) + y(n-3) + ... + y2 + y1
-----------------------------------------------
=  N.Y(n) + (N-1).Y(n-1) + (N-2).Y(n-2) + (N-3).Y(n-3) + ... + 2.y2 + 1.Y1

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
int Sum = 0, Sxy = 0;
 
for i=n downto 1
  Sum += y[i]
  Sxy += Sum
next

**progfou** · 15/02/2011, 23h24

pseudocode ta méthode est top !

Ca ne fait que des additions, ça me plait bien.
Reste une seule chose à régler : les outliers.
En fait, je sais que ça ne varie pas énormément, donc une technique type passage en revue du tableau et détection de l'outlier par sa variance correspond mais ajoute un peu de complexité (un tour sup du tableau).
Il faut que je vois ce qui est le plus efficace en comparaison d'un filtre y(n) = (y(n-1)+y(n))/2

**pseudocode** · 16/02/2011, 00h17

Envoyé par progfou

pseudocode ta méthode est top !

Ca ne fait que des additions, ça me plait bien.
Reste une seule chose à régler : les outliers.
En fait, je sais que ça ne varie pas énormément, donc une technique type passage en revue du tableau et détection de l'outlier par sa variance correspond mais ajoute un peu de complexité (un tour sup du tableau).
Il faut que je vois ce qui est le plus efficace en comparaison d'un filtre y(n) = (y(n-1)+y(n))/2

utiliser le médian plutôt que la moyenne, ca serait mieux.

Code java :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
public static void main(String[] args) {
 
	int[] y = {31134,31288,31175,31318,31177,31368,31240,31363,31278,31434,31284,31463,31375,31511,31366,31527,31411,31589,31444,31601,31501,31653,31523,31702,31574,31734,31594,31767,31643,31814};
 
	int N = y.length-2;
	int Sx = N*(N+1)/2;
	int Sxx = N*(N+1)*(2*N+1)/6;
	int Sy = 0, Sxy = 0;
 
	for (int i = N; i >= 1; i--) {
		Sy += median(y[i-1],y[i],y[i+1]);
		Sxy += Sy;
	}
 
	double a = (double)(N*Sxy - Sx*Sy) / (N*Sxx - Sx*Sx);
	double b = (Sy-a*Sx)/N;
	System.out.printf("y = %.2f.x + %.2f",a,b);
}
 
static int median(int a,int b,int c) {
	if (a<b && b<c) return b;
	if (a>b && b>c) return b;
	if (b<a && a<c) return a;
	if (b>a && a>c) return a;
	return c;
}