Régression linéaire par morceaux

Version imprimable

23/06/2017, 18h59
Négrier

Régression linéaire par morceaux

Connaitriez vous un algorithme ou un programme (pascal ou VB) de régression linéaire par morceaux?
23/06/2017, 19h22
souviron34

euh..

Peux-tu expliquer un peu ?

Parce que en général ça s'appelle juste "interpolation linéaire" : entre chaque couple de points on interpole linéairement..

Donc le programmer est extrêmement simple.... dans n'importe quel langage...

On a une base de points initiaux (un tableau) Pts (x,y)

Pour un nouveau point, on cherche entre quels indices il se situe. (i, i+1)

Puis on fait une règle de trois.. On "interpole linéairement"
23/06/2017, 20h09
dourouc05

Donc, on peut faire une régression plus complexe que linéaire (quadratique, polynomiale en général), puis en faire une approximation linéaire par morceaux.
24/06/2017, 11h30
pseudocode

Je pense que sa question c'est "comment déterminer les morceaux" pour avoir une approximation optimale. :)

Si c'est bien cela, il y a plusieurs méthodes dans la littérature (optimal piecewise linear approximation).

Le plus simple (mais pas forcément optimal) c'est de faire une suite d'approximations locales:
1. Approximer la courbe par un segment qui rejoint les 2 extrémités de la courbe: [A,B].
2. Chercher le point P de la courbe qui est le plus éloigné de l'approximation (= le plus loin du segment).
3. Remplacer le segment [A,B] par les deux segments [A, P] et [P, B]
4. Recommencer à l'étape 2 jusqu'à atteindre un erreur acceptable.

Sinon il y a la méthode inverse dans laquelle on fait un découpage uniforme de la courbe en petits segments,
puis on on fusionne itérativement les segments jusqu'à atteindre un erreur acceptable (Douglas-Peucker algorithm)
26/06/2017, 00h05
Vincent PETIT

Salut,

Citation:

Envoyé par souviron34

Parce que en général ça s'appelle juste "interpolation linéaire"

Pour moi une régression linéaire c'est chercher une droite qui traversent au mieux l'ensemble des points alors qu'une interpolation linéaire c'est chercher une courbe qui passe par tous les points.

Non ?
26/06/2017, 05h56
souviron34

Si, je suis d'accord....

MAIS

le problème est le "par morceaux"

J'ai bien peur que le PO n'aie pas très très bien compris.. -ou expliqué...

;)
26/06/2017, 09h34
wiwaxia

Régression linéaire par morceaux

Salut, :D

pseudocode a apporté la réponse la plus claire.

La seule information supplémentaire que l'on pourrait souhaiter, c'est l'aspect général du nuage de points: la "courbe moyenne" présente-t-elle une courbure continue, ou des parties quasi-rectilignes ? Les domaines de régression linéaire sont-ils éventuellement prédéfinis en limite ou en nombre ? Cela pourrait orienter la recherche du meilleur algorithme.
26/06/2017, 11h30
Négrier

1 pièce(s) jointe(s)

Détail sur le problème précis (jeu de données)

j'ai une fonction Y=f(X) définie sur un très grand nombre de points (27000) lorsque je représente cette fonction en
Excel j'obtiens une ligne brisée de quelques segments (entre 5 et 10). Mais ça c'est ce que je constate
sur mon graphe mais je ne connais pas d'algorithme/fonction qui me permettrait de récupérer, à partir de tous ces points, la suite
des segments décrits par leurs extrémités (Xdébut, Ydébut) --> (Xfin, Yfin).

Je suis à la recherche code ou simplement un algorithme qui me permettrait de faire cela ?
26/06/2017, 12h48
wiwaxia

1 pièce(s) jointe(s)

Régression linéaire par morceaux

Alors la solution me paraît toute trouvée: le calcul de la pente locale sur les (2*n) points encadrant le point considéré, de rang (k), par la formule:

p_k = (n+1)^-1*[(y_k+1 + y_k+2 + ... + y_k+n) - (y_k-1 + y_k-2 + ... + y_k-n)] sur le domaine [n ; N - k].

Prendre une valeur modérée pour (n): 3 à 10 ? à toi de voir; il n'est même pas nécessaire de tenir compte du facteur constant (n+1)^-1; la grandeur calculée sera alors p'_k = (n+1) * p_k .

L'essentiel est que le graphe de la suite y_(k) = p_k présente une série de paliers horizontaux séparés par une variation très rapide au voisinage des points de raccordement.

Pièce jointe 289078

Les bornes apparaîtront d'autant plus nettement que (n) est plus petit ... d'où l'intérêt, pour ce dernier, de ne pas lui donner une valeur trop élevée - ce qui favorisera aussi la rapidité du calcul.
26/06/2017, 14h44
wiwaxia

Régression linéaire par morceaux

Autre solution, encore plus directe: les points de raccordement des parties quasi-rectilignes (qu'il ne sera peut-être pas nécessaire d'incorporer aux calculs de régression linéaire) se distingueront par un extremum de la variation de la pente calculée sur (2*n) points deux à deux symétriquement disposés, par une combinaison linéaire semblable à la précédente:

q_k = (p_k+1 + p_k+2 + ... + p_k+n) - (p_k-1 + p_k-2 + ... + p_k-n)

On obtient ainsi pour (n = 3): q_k = S_i=-6⁶(C_i * y_k+i)
avec: C₀ = -6 , C₁ = C_-1 = -4 , C₂ = C_-2 = -1 , C₃ = C_-3 = 2 ,
C₄ = C_-4 = 3 , C₅ = C_-5 = 2 , C₆ = C_-6 = 1 ;
le pic observé, centré en (k), présente une demi-largeur à la base égale à 3.
26/06/2017, 17h05
dourouc05

Sinon, fixer les bornes à la main, ne considérer que les points de ce que tu veux comme droite, faire une régression linéaire sur ce sous-ensemble de points.
26/06/2017, 19h51
Négrier

récupérer une information synthétique

mon objectif est de pouvoir synthétiser l'information effectivement utile de mes 27000 points en une dizaine de segments.
L'objet informatique ainsi obtenu est beaucoup plus facile à manipuler et il contient l'essentiel de l'information.
26/06/2017, 23h09
tbc92

Sur l'exemple fourni, les données sont presque parfaitement alignées. Si tes données sont systématiquement aussi "propres", la solution de Wiwaxia va fonctionner.
Etape 1 : Remplacer chaque nombre par une moyenne mobile. (une moyenne mobile, c'est la moyenne par exemple du nombre en cours, des 10 nombres précédents et des 10 nombres suivants, c.a.dire somme des 21 nombres, divisée par 21) . On fait l'impasse sur les 10 premiers nombres , et sur les 10 derniers. On a donc une série de 27000-20 = 26980 moyennes mobiles MM.

Etape 2 : Calculer les écarts EC(n) = MM(n) - MM(n-1) ; ici, on voit nettement un premier groupe où les écarts EC sont très proches de -0.852 , puis un 2nd groupe où les écarts sont proches de -0.48 etc etc.. Tant que EC(n) est proche de EC(n-1) avec moins de 1% d'écart (seuil qu'on peut affiner), on considère qu'on est sur la même portion de droite.
Pour notre premier segment , on va prendre par exemple un point vers le début (-14.825, 16 ) et un point vers la fin ( -1488.217, 1748) Et on trouve l'équation de la droite qui passe par ces 2 points : Y = -0.8506822 x + 1.213988
On peut aussi prendre les moyennes mobiles, au lieu des points d'origine, pour limiter le risque de tomber sur un point atypique.

Ensuite, on a une zone où notre écart EC n'est pas stable (en fait il n'est pas stable parce qu'on a pris des moyenne mobiles ; si on travaillait avec les données d'origine, l'Ecart EC serait plus stable ici). Puis on va bien identifier un long segment avec une pente de -0.48. Idem, on prend un point vers le début de ce segment et un point vers la fin, et on sait calculer l'équation de cette droite.

Et pour la zone intermédiaire, on calcule où nos 2 droites se coupent. Avec des données très propres comme celles que tu as, ça ne devrait pas poser de problème.
26/06/2017, 23h53
wiwaxia

Régression linéaire par morceaux
Afin d'éviter tout imbroglio, je signale avoir modifié (et simplifié) le contenu calculatoire de mon précédent message (#10).

Les variations au premier et second ordre s'expriment finalement par la même formule, utilisant n = 3 paires de termes.

Voici ce que l'on obtient au voisinage d'une discontinuité définie par: u_(k<0) = 0; u₍₀₎ = 5 ; u_(k>0) = 10 .
Code:

1 2 3 k-4 k-3 k-3 k-1 k k+1 k+2 k+3 k+4 0 5 15 25 30 25 15 5 0
Chaque rupture de pente est ainsi repérable par un pic de hauteur proportionnelle à la discontinuité étudiée.

Douglas Peucker

Merci pour toutes vos réponses.

Le message de pseudocode m'a orienté vers l'algorithme de Douglas Peucker (je me doutais qu'il devait y avoir des références sur
cette "régression linéaire par morceaux" mais je ne savais pas par où commencer). Je vais donc regarder l'algo ci-dessous en plus de tester
les solutions que vous m'avez proposées.

Code:

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
27
28
29
function DouglasPeucker(PointList[], epsilon)
  // Trouve le point le plus éloigné du segment
  dmax = 0
  index = 0
  for i = 2 to (length(PointList) - 1)
    d = DistancePointSegment(PointList[i], Segment(PointList[1], PointList[end])) 
    if d > dmax
      index = i
      dmax = d
    end
  end
 
  // Si la distance dmax est supérieure au seuil, on simplifie
  if dmax > epsilon
    // Appel récursif de la fonction
    recResults1[] = DouglasPeucker(PointList[1…index], epsilon)
    recResults2[] = DouglasPeucker(PointList[index…end], epsilon)
 
    // Construit la liste des résultats à partir des résultats partiels
    ResultList[] = {recResults1[1…end-1] recResults2[1…end]}
 
  else
    // Tous les points sont proches → renvoie un segment avec les extrémités
    ResultList[] = {PointList[1], PointList[end]}
  end
 
  // Renvoie le nouveau résultat
  return ResultList[]
end

4 pièce(s) jointe(s)

L’algorithme de Douglas-Peuckertraduit en langage Pascal

L’algorithme de Douglas-Peucker donne de très bons résultats, la version itérative du programme associé est très efficace.
Le choix du bon seuil (« epsilonLocal » dans le code ci-dessous) est bien entendu crucial.

C'est en gros la méthode suggérée par pseudocode

Exemple 1, résultat obtenu : 27687 points (graphe de gauche) synthétisé en à 6 points (graphe de droite)

Pièce jointe 289853 Pièce jointe 289856

Exemple 2, résultat obtenu : 13128 points (graphe de gauche) synthétisé en à 7 points (graphe de droite)

Pièce jointe 289861 Pièce jointe 289863

Le codage en pascal de l’algorithme de Douglas-Peucker dans sa version itérative (Ramer)

Code:

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
 
 
 
{
 
Algorithme de  Ramer-Douglas-Peucker
À chaque étape on parcourt tous les nœuds entre les bornes et on sélectionne
le nœud le plus éloigné du segment formé par les bornes :
 
1.     s'il n'y a aucun nœud entre les bornes l'algorithme se termine,
2.     si cette distance est inférieure à un certain seuil on supprime tous les nœuds entre les bornes,
3.     si elle est supérieure la polyligne n'est pas directement simplifiable.
  On appelle de manière récursive l'algorithme sur deux sous-parties de la
  polyligne : de la première borne au nœud distant, et du nœud distant à la borne finale.
 
 
Iterative version of Ramer-Douglas-Peucker line-simplification algorithm
 
}
 
 
 
unit uAlgoDouglasPucker;
 
interface
 
uses system.generics.collections, system.types, system.Math.vectors,
  system.classes;
 
type
 
  TListePoints = array of TpointF;
 
 
procedure DouglasPeucker(const points: TListePoints; epsilonLocal: double;
  var res: TListePoints);
 
 
implementation
 
type
 
  TKeyValuePair = record
    key, value: integer;
    constructor create(aun, adeux: integer);
  end;
 
 
 
function DouglasPeuckerInterne(const points: TListePoints;
  startIndex, lastIndex: integer; epsilonLocal: double)
  : TBooleanDynArray; forward;
 
function pointLineDistance(point, start, fin: TpointF): double; forward;
 
 
 
procedure DouglasPeucker(const points: TListePoints; epsilonLocal: double;
  var res: TListePoints);
 
var
  bitArray: TBooleanDynArray;
  i, hg: integer;
 
begin
  bitArray := DouglasPeuckerInterne(points, 0, high(points), epsilonLocal);
  res := nil;
  hg := high(res);
  for i := Low(points) to High(points) do
  begin
    if bitArray[i] then
    begin
      hg := hg + 1;
      setLength(res, hg + 1);
      res[hg] := points[i];
    end;
  end;
end;
 
 
 
function DouglasPeuckerInterne(const points: TListePoints;
  startIndex, lastIndex: integer; epsilonLocal: double): TBooleanDynArray;
var
  pile: TStack<TKeyValuePair>;
  bitArray: TBooleanDynArray;
  dmax, d: double;
  index, i, globalStartIndex: integer;
begin
  globalStartIndex := startIndex;
  setLength(bitArray, lastIndex - startIndex + 1);
  for i := Low(bitArray) to High(bitArray) do
    bitArray[i] := TRUE;
 
  pile := TStack<TKeyValuePair>.create();
  pile.Push(TKeyValuePair.create(startIndex, lastIndex));
  while (pile.Count > 0) do
  begin
    startIndex := pile.Peek.key;
    lastIndex := pile.Peek.value;
    pile.Pop;
    dmax := 0.0;
    index := startIndex;
    for i := index + 1 to lastIndex - 1 do
    begin
      if bitArray[i - globalStartIndex] then
      begin
        d := pointLineDistance(points[i], points[startIndex],
          points[lastIndex]);
        if d > dmax then
        begin
          index := i;
          dmax := d;
        end;
      end;
    end;
 
    if dmax > epsilonLocal then
    begin
      pile.Push(TKeyValuePair.create(startIndex, index));
      pile.Push(TKeyValuePair.create(Index, lastIndex));
    end
    else
    begin
      for i := startIndex + 1 to lastIndex - 1 do
        bitArray[i - globalStartIndex] := FALSE;
    end;
  end;
 
  result := bitArray;
 
end;
 
 
 
function pointLineDistance(point, start, fin: TpointF): double;
var
  n, d: double;
begin
  if (start.X = fin.X) and (start.Y = fin.Y) then
    result := point.Distance(start)
  else
  begin
    n := Abs((fin.X - start.X) * (start.Y - point.Y) - (start.X - point.X) *
      (fin.Y - start.Y));
    d := Sqrt((fin.X - start.X) * (fin.X - start.X) + (fin.Y - start.Y) *
      (fin.Y - start.Y));
    result := n / d;
  end;
end;
 
 
 
{ TKeyValuePair }
 
constructor TKeyValuePair.create(aun, adeux: integer);
begin
  key := aun;
  value := adeux;
end;
 
end.