CURVE FITTER WANTED!
Author Message
CURVE FITTER WANTED!

am looking for a Pascal program (source if possible) which does this:

I have a polynom ( a*x^4+b*x^3+c*x^2+d*x+e) and x,y data and it determines
a,b,c,d,e so the curve fits the x,y data using less square error method.

I need it for 2-5 degree polynoms.
( from a*x^2+b*x+c, to a*x^5+b*x^4+c*x^3+d*x^2+e*x+f )

Thanks
Joe

Wed, 18 Jun 1902 08:00:00 GMT
CURVE FITTER WANTED!

Quote:

>  am looking for a PASCAL program (source if possible) which does this:

> I have a polynom ( a*x^4+b*x^3+c*x^2+d*x+e) and x,y data and it determines
> a,b,c,d,e so the curve fits the x,y data using less square error method.

<snip>

Quote:
> Thanks
> Joe

Timo Salmi just posted the following link:

303778 May 2 1991 ftp://garbo.uwasa.fi/pc/turbopas/nrpas13.zip
nrpas13.zip Numerical Recipes Pascal shareware version, 303K!

The Numerical Recipes should have this for you.
Also, see the book, _Numerical Recipes in Pascal_ by
Press, Flannery, et al.

Al Moore

Wed, 18 Jun 1902 08:00:00 GMT
CURVE FITTER WANTED!

:I have a polynom ( a*x^4+b*x^3+c*x^2+d*x+e) and x,y data and it determines
:a,b,c,d,e so the curve fits the x,y data using less square error method.

Just in case it isn't a typo. "The least squares method".

All the best, Timo

....................................................................

Moderating at ftp:// & http://garbo.uwasa.fi/ archives 193.166.120.5
Department of Accounting and Business Finance  ; University of Vaasa

Spam foiling in effect.  My email filter autoresponder will return a
required email password to users not yet in the privileges database.
Advice on spam foiling at http://www.uwasa.fi/~ts/info/spamfoil.html

Wed, 18 Jun 1902 08:00:00 GMT
CURVE FITTER WANTED!

Quote:

> am looking for a PASCAL program (source if possible) which does this:

>I have a polynom ( a*x^4+b*x^3+c*x^2+d*x+e) and x,y data and it determines
>a,b,c,d,e so the curve fits the x,y data using less square error method.

Try
www.efg2.com/lab/library/Delphi/MathFunctions/StatisticsAndProbabilit...

efg
_________________________________
efg's Computer Lab:       www.efg2.com/lab
Delphi Books:  www.efg2.com/lab/TechBooks/Delphi.htm

Overland Park, KS  USA

Wed, 18 Jun 1902 08:00:00 GMT
CURVE FITTER WANTED!

Quote:

> am looking for a PASCAL program (source if possible) which does this:

>I have a polynom ( a*x^4+b*x^3+c*x^2+d*x+e) and x,y data and it determines
>a,b,c,d,e so the curve fits the x,y data using less square error method.

>I need it for 2-5 degree polynoms.
>( from a*x^2+b*x+c, to a*x^5+b*x^4+c*x^3+d*x^2+e*x+f )

>Thanks

Here is some code I produced some years ago.
Without some mathematial knowledge it's hard to
understand but you can use it as a 'black box'.

Best Regards
R.Fischer
(hope the posting is not too long for usenet)

----------------------------------------------
PROGRAM LineFit;
Uses    Crt, Graph;       { attention with Crt on fast machines }

CONST   MaxK   = 25;      { max. number of fitting functions }
MaxN   = 2000;    { max. number of xy-points }

TYPE    Real = System.Double;   { it's faster and more exact }
Vector = ARRAY[1..MaxN] OF Real;
Matrix = ARRAY[1..MaxK,1..MaxK] OF Real;

PROCEDURE Gaussj(VAR a: matrix; n: integer; VAR b: vector);
{- Solving a LinEqu. (Gauss-Jordan from Numerical Recipes for Pascal) -}
TYPE  glnp = ARRAY [1..MaxK] OF integer;
VAR
big,dum,pivinv: real;
i,icol,irow,j,k,l,ll: integer;
indxc,indxr,ipiv: glnp;
BEGIN
FOR j := 1 TO n DO BEGIN
ipiv[j] := 0
END;
FOR i := 1 TO n DO BEGIN
big := 0.0;
FOR j := 1 TO n DO BEGIN
IF (ipiv[j] <> 1) THEN BEGIN
FOR k := 1 TO n DO BEGIN
IF (ipiv[k] = 0) THEN BEGIN
IF (abs(a[j,k]) >= big) THEN BEGIN
big := abs(a[j,k]);
irow := j;
icol := k
END
END ELSE IF (ipiv[k] > 1) THEN BEGIN
writeln('pause 1 in GAUSSJ - singular matrix'); readln
END
END
END
END;
ipiv[icol] := ipiv[icol]+1;
IF (irow <> icol) THEN BEGIN
FOR l := 1 TO n DO BEGIN
dum := a[irow,l];
a[irow,l] := a[icol,l];
a[icol,l] := dum
END;
FOR l := 1 TO 1 DO BEGIN
dum := b[irow];
b[irow] := b[icol];
b[icol] := dum
END
END;
indxr[i] := irow;
indxc[i] := icol;
IF (a[icol,icol] = 0.0) THEN BEGIN
writeln('pause 2 in GAUSSJ - singular matrix'); readln
END;
pivinv := 1.0/a[icol,icol];
a[icol,icol] := 1.0;
FOR l := 1 TO n DO BEGIN
a[icol,l] := a[icol,l]*pivinv
END;
FOR l := 1 TO 1 DO BEGIN
b[icol] := b[icol]*pivinv
END;
FOR ll := 1 TO n DO BEGIN
IF (ll <> icol) THEN BEGIN
dum := a[ll,icol];
a[ll,icol] := 0.0;
FOR l := 1 TO n DO BEGIN
a[ll,l] := a[ll,l]-a[icol,l]*dum
END;
FOR l := 1 TO 1 DO BEGIN
b[ll] := b[ll]-b[icol]*dum
END
END
END
END;
FOR l := n DOWNTO 1 DO BEGIN
IF (indxr[l] <> indxc[l]) THEN BEGIN
FOR k := 1 TO n DO BEGIN
dum := a[k,indxr[l]];
a[k,indxr[l]] := a[k,indxc[l]];
a[k,indxc[l]] := dum
END
END
END
END;   {--- Gaussj ---}

PROCEDURE FCT(k:Integer; xarg: Real; VAR PHI: Vector);
{ the set of the fitting functions }
{ for example Polynomials }
VAR   i  : Integer;
BEGIN
PHI[1]:=1;    { polynomials }
FOR i:=2 TO k DO PHI[i]:=PHI[i-1]*xarg;

{ Another set of 'better' (for the fitting problem) functions can be used too }
{ These functions must be 'linearly independent' (polynomials are sure) }
{ for example:
PHI[1]:=1;
PHI[2]:=1/(xarg+1);
PHI[3]:=xarg/(xarg*xarg+1);
PHI[4]:=sin(xarg);
PHI[5]:=exp(xarg);
....
....
}
END;   {---- FCT ----}

FUNCTION CalcFitting(k : Integer; A: Vector; xarg : Real) : Real;
{ Calc one Y-Point of the fitting curve }
VAR   PHI  : Vector;
sum  : Real;
i    : Integer;
BEGIN
FCT(k,xarg,PHI);
sum:=0;
FOR i:=1 TO k DO sum:=sum + A[i]*PHI[i];
CalcFitting:=sum;
END;

PROCEDURE DoFitting(n,k : Integer; X,Y : Vector; VAR A : Vector);
{ Calc the coefficients of the fitting system - exported in A}
TYPE VectPoint   = ^Vector;
FType       = ARRAY[1..MaxK] OF VectPoint;
VAR  i, j, l     : Integer;
PHI         : Vector;
F           : FType;
M           : ^Matrix;
BEGIN
New(M);                      { Reserve memory on the heap }
FOR i:=1 TO k DO New(F[i]);  { Missing: is enough there ? }

FOR i:=1 TO n DO         { Prepare the approximation }
BEGIN
FCT(k,X[i],PHI);
FOR j:=1 TO k DO F[j]^[i]:=PHI[j];
END;
FOR j:=1 TO k DO         { Build the matrix of the LinEqu. }
FOR i:=1 TO k DO
BEGIN  M^[i,j]:=0;
FOR l:=1 TO N DO
M^[i,j]:=M^[i,j] +  F[i]^[l]*F[j]^[l];
M^[j,i]:=M^[i,j];
END;
FOR i:=1 TO k DO         { Build the right side of the LinEqu.}
BEGIN  A[i]:=0;
FOR l:=1 TO n DO A[i]:=A[i] + F[i]^[l]*Y[l];
END;

Gaussj(M^,k,A);    { Solve the system }
{
Missing:
IF (Matrix is singular) THEN Errorhandling;
}
FOR i:=1 TO k DO Dispose(F[i]);   { Clear reserverd mem }
Dispose(M);
END;   {---- DoFitting ----}

VAR   X, Y, A                          : Vector;
n, i, k, GraphDriver, GraphMode  : Integer;
delta, Zero, xpos, ypos          : Real;
ch                               : Char;
st                               : String[3];

BEGIN     { ---- MAIN   (example)  ---- }
n:=300;      { 300 data-points }
GraphDriver:=Detect;
InitGraph(GraphDriver,GraphMode,'');

delta:=(GetMaxX)/(n-1);  Zero:=GetMaxY/2;
FOR i:=1 TO n DO      { Build an example }
BEGIN
x[i]:=(i-1)*delta;
y[i]:=Zero + (Sin(X[i]/50)*Zero*0.5)*(X[i]/400) + (Random-0.5)*25 ;
END;

delta:=GetMaxX/500;    { Show 15 line-fittings with polynomials }
FOR k:=1 TO 15 DO      { k = (degree of the polynomial)-1 }
BEGIN
ClearDevice;
Str(k-1:1,st);  SetColor(White);
OutTextXY(10,10,'Degree of the fitting polynomial: '+st);
SetColor(LightRed);          { Show the Points }
FOR i:=1 TO n DO Circle(Round(x[i]),Round(y[i]),2);
DoFitting(n,k, X,Y, A); { Do the calculation of the fitting polynomial }
FOR i:=0 TO 500 do      { Show the fitting curve with 501 points }
BEGIN
xpos:=i*delta;
ypos:=CalcFitting(k,A,xpos);
PutPixel(Round(xpos),Round(ypos),White);
END;
END;
CloseGraph;
END.

Wed, 18 Jun 1902 08:00:00 GMT

 Page 1 of 1 [ 5 post ]

Relevant Pages