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:

unit GaussJordan_Unit;

interface


Type TTMatrix = Array of Array of Single;
Type TTVector = Array of Single;
Function GaussJordan(mat:TTMatrix; var x:TTVector):Boolean;

implementation


Function GaussJordan(mat:TTMatrix; var x:TTVector):Boolean;
var
     n         : Integer;
     i,j,k     : Integer;
     maxrow    : Integer;
     tmp       : Single;
begin

     n := Length(mat);

     for i:= 0 to n-1 do
     begin
          //Find the row with the largest first value
          maxrow := i;
          for j := i+1 to n-1 do
          begin
               if abs(mat[j,i]) > abs(mat[maxrow,i]) then
                    maxrow := j;
          end;

          //Swap the maxrow and ith row
          for k := i to n do
          begin
               tmp := mat[i,k];
               mat[i,k] := mat[maxrow,k];
               mat[maxrow,k] := tmp;
          end;

          //Singular matrix?
          if abs(mat[i][i]) < 0.0000001 then
          begin
               Result := False;
               Exit;
          end;

          //Eliminate the ith element of the jth row
          for j := i+1 to n-1 do
          begin
               for k := n downto i do
               begin
                    mat[j,k] := mat[j,k] - (mat[i,k] * mat[j,i] / mat[i,i]);
               end;
          end;
     end;

     //Do the back substitution
     setLength(x, n);
     for j := n-1 downto 0 do
     begin
          tmp := 0;

          for k := j+1 to n-1 do
               tmp := tmp + mat[j,k] * x[k];

          x[j] := (mat[j,n] - tmp) / mat[j,j];
     end;

     Result := True;
end;

end.					

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:

unit LinReg_Unit;

interface

Type TTLRResultsG1 = packed record
     a,b       : Single;
end;

Type TTLRResultsG2 = packed record
     a,b,c     : Single;
end;

Type TTLRPoint = packed record
     x,y       : Single;
end;

Type TTLRInputs = packed record
     data      : Array of TTLRPoint;
     count     : Cardinal;
end;

Function DoLinRegG1(points:TTLRInputs; out r:TTLRResultsG1):Boolean;
Function DoLinRegG2(points:TTLRInputs; out r:TTLRResultsG2):Boolean;

implementation

uses
     GaussJordan_Unit;


Function DoLinRegG1(points:TTLRInputs; out r:TTLRResultsG1):Boolean;
var
     mat       : TTMatrix;
     x         : TTVector;
     i         : Integer;
begin
     SetLength(mat, 2, 3);

     //Für Regression einer Parabel: f(x)=a*x + b
     //[x^2]  [x]  [y·x]
     //[x]  [1]  [y]

     //Gaußsummen [x^2] = p1.x^2 + p2.x^2 + ... + pn.x^2

     mat[0,0] := 0; mat[0,1] := 0; mat[0,2] := 0;
     mat[1,0] := 0; mat[1,1] := 0; mat[1,2] := 0;

     for i := 0 to points.count-1 do
     begin
          mat[0,0] := mat[0,0] + points.data[i].x*points.data[i].x;

          mat[0,1] := mat[0,1] + points.data[i].x;
          mat[1,0] := mat[0,1];

          mat[1,1] := mat[1,1] + 1;

          mat[0,2] := mat[0,2] + points.data[i].y*points.data[i].x;
          mat[1,2] := mat[1,2] + points.data[i].y;
     end;

     if not GaussJordan(mat, x) then
     begin
          Result := False;
          Exit;
     end;

     r.a := x[0];
     r.b := x[1];

     Result := True;
end;

Function DoLinRegG2(points:TTLRInputs; out r:TTLRResultsG2):Boolean;
var
     mat       : TTMatrix;
     x         : TTVector;
     i         : Integer;
begin
     SetLength(mat, 3, 4);

     //Für Regression einer Parabel: f(x)=a*x^2 + b*x + c
     //[x^4]  [x^3]  [x^2]     [y·x^2]
     //[x^3]  [x^2]  [x]       [y·x]
     //[x^2]  [x]       [1]       [y]

     //Gaußsummen [x^2] = p1.x^2 + p2.x^2 + ... + pn.x^2

     mat[0,0] := 0; mat[0,1] := 0; mat[0,2] := 0; mat[0,3] := 0;
     mat[1,0] := 0; mat[1,1] := 0; mat[1,2] := 0; mat[1,3] := 0;
     mat[2,0] := 0; mat[2,1] := 0; mat[2,2] := 0; mat[2,3] := 0;

     for i := 0 to points.count-1 do
     begin
          mat[0,0] := mat[0,0] + points.data[i].x*points.data[i].x*
                                 points.data[i].x*points.data[i].x;

          mat[0,1] := mat[0,1] + points.data[i].x*points.data[i].x*points.data[i].x;
          mat[1,0] := mat[0,1];

          mat[0,2] := mat[0,2] + points.data[i].x*points.data[i].x;
          mat[1,1] := mat[0,2];
          mat[2,0] := mat[0,2];

          mat[1,2] := mat[1,2] + points.data[i].x;
          mat[2,1] := mat[1,2];

          mat[2,2] := mat[2,2] + 1;

          mat[0,3] := mat[0,3] + points.data[i].y*(points.data[i].x*points.data[i].x);
          mat[1,3] := mat[1,3] + points.data[i].y*points.data[i].x;
          mat[2,3] := mat[2,3] + points.data[i].y;
     end;

     if not GaussJordan(mat, x) then
     begin
          Result := False;
          Exit;
     end;

     r.a := x[0];
     r.b := x[1];
     r.c := x[2];

     Result := True;
end;

end.					

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:

procedure TForm1.FormCreate(Sender: TObject);
var
     mat       : TTMatrix;
     x         : TTVector;
     points    : TTLRInputs;
     r2        : TTLRResultsG2;
     r1        : TTLRResultsG1;
begin
     //SetLength(mat, 3, 4);
     //mat[0,0] := 1; mat[0,1] := 2; mat[0,2] := 1; mat[0,3] := 1;
     //mat[1,0] := 4; mat[1,1] := 5; mat[1,2] := 2; mat[1,3] := 2;
     //mat[2,0] := 1; mat[2,1] := 2; mat[2,2] := 3; mat[2,3] := 2;
     //GaussJordan(mat, x);
     //ShowMessage(FloatToStr(x[0]) + ' ' +  FloatToStr(x[1]) + ' ' + FloatToStr(x[2]));

     SetLength(points.data, 4);
     points.count := 4;

     points.data[0].x := 2;
     points.data[0].y := -3;

     points.data[1].x := 3;
     points.data[1].y := 0;

     points.data[2].x := 6;
     points.data[2].y := 7;

     points.data[3].x := -1;
     points.data[3].y := -4;

     DoLinRegG1(points, r1);
     DoLinRegG2(points, r2);

     ShowMessage('f(x) = ' + FloatToStr(r1.a) + '*x + ' +  FloatToStr(r1.b));
     ShowMessage('f(x) = ' + FloatToStr(r2.a) + 'x^2 + ' +  FloatToStr(r2.b) + '*x + ' + FloatToStr(r2.c));
end;