{$mode objfpc}{$m+}{$N-}
program Quadratic_fit;
uses
Classes,
Math;
var
i, j, k: integer;
x: array of double = (27.23103448, 27.22068966,
27.25862069, 27.23793103, 27.12068966, 27.21034483, 27.59310345,
28.17931034, 28.54137931, 28.95517241, 29.45517241, 30.63793103,
30.17931034, 30.56551724, 30.13103448, 29.95172414, 29.68275862,
29.06206897, 28.38275862, 28.03793103, 27.66551724, 27.53448276, 27.43793103);
y: array of double = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18, 19,20,21,22,23);
coeffs: array of double;
type
TLocalVector = array of double;//real;
TPolyReg = class //PolynomialRegression
private
n, np1, np2, tnp1: longint;
N1: longint;
x1, y1, A: TLocalVector;
B: array of array of double;
tmp, t: double;
Len_x, Len_y, Len_a: longint;
nm1: integer;
Result: boolean;
public
constructor init(); overload;
destructor done(); virtual;
function fitIt(const x: TLocalVector; const y: TLocalVector;
const order: integer; var coeffs: TLocalVector): boolean;
end;
constructor TPolyReg.init;
begin
inherited Create;
end;
destructor TPolyReg.done;
begin
inherited Destroy;
end;
function TPolyReg.fitIt(const x: TLocalVector; const y: TLocalVector;
const order: integer; var coeffs: TLocalVector): boolean;
begin
Result := True;
// The size of xValues and yValues should be same. (
if Length(x) <> Length(y) then
begin
writeln('Kich co cua mang X va mang Y khong bang nhau');
Result := False;
end;
// The size of xValues and yValues cannot be 0, should not happen
Len_x := Length(x);
Len_y := Length(y);
if (Len_x = 0) or (Len_y = 0) then
begin
writeln('Kich co cua mang X hoac mang Y bang 0');
Result := False;
end;
n := order;
np1 := n + 1;
np2 := n + 2;
tnp1 := 2 * n + 1;
// X = vector that stores values of sigma(xi^2n).
SetLength(x1, tnp1);
for i := 0 to tnp1 - 1 do
begin
x1[i] := 0;
for j := 0 to Len_x - 1 do
x1[i] += power(x[j], i);
end;
// a = vector to store final coefficients.
SetLength(a, np1);
// B = normal augmented matrix that stores the equations.
SetLength(B, np1, np2); // n x m = np1 x np2
for i := 0 to n do
for j := 0 to n do
B[i][j] := x1[i + j];
// Y = vector to store values of sigma(xi^n * yi).
SetLength(Y1, np1);
for i := 0 to np1 - 1 do
begin
y1[i] := 0;
for j := 0 to n1 - 1 do
y1[i] += power(x[j], i) * y[i];
end;
// Load values of Y as last column of B.
for i := 0 to n do
B[i][np1] := y1[i];
n += 1;
nm1 := n - 1;
// Pivotisation of the B matrix.
for i := 0 to n - 1 do
for k := i + 1 to n - 1 do
begin
if (B[i][i] < B[k][i]) then
for j := 0 to n do
begin
tmp := B[i][j];
B[i][j] := B[k][j];
B[k][j] := tmp;
end;
end;
// Performs the Gaussian elimination.
// (1) Make all elements below the pivot equals to zero
// or eliminate the variable.
for i := 0 to nm1 - 1 do
for k := 0 to n - 1 do
begin
t := B[k][i] / B[i][i];
for j := 0 to n do
B[k][i] -= t * B[i][j]; // (1)
end;
// Back substitution.
// (1) Set the variable as the rhs of last equation
// (2) Subtract all lhs values except the target coefficient.
// (3) Divide rhs by coefficient of variable being calculated.
{I am not sure about this correction
n:=n-1;//engkin
nm1:=n-1;//engkin}
for i := nm1 downto 0 do
begin
a[i] := B[i][n]; // (1)
for j := 0 to n - 1 do
if j <> i then
a[i] -= B[i][j] * a[j]; // (2)
a[i] /= B[i][i]; // (3)
end;
Len_a := Length(a);
SetLength(coeffs, Len_a);
for i := 0 to Length(a) - 1 do
coeffs[i] := a[i];
end;
var
PolyReg: TPolyReg;
begin
PolyReg := TPolyReg.init;
WriteLn('fitIt: ',PolyReg.fitIt(x, y, 10, coeffs));
WriteLn('coefficients:');
for j := 0 to High(coeffs) do
WriteLn(' ', coeffs[j]);
writeln;
PolyReg.done;
ReadLn;
end.
Recent
