unit LeastSquares;

{Public-domain code for solving least-squares problems.}

{ Compatible with IMAGE (an image-measurement program by Wayne Rasband at NIH)}
{ and a functional substitute for the TurboPascalª Numerical Methods Toolbox in}
{ that application.  Not all calling-options of the Numerical Methods package}
{ are implemented.}

{ Written September 1989 under MPW Pascal by}
{    Dr. John C. Russ}
{    North Carolina State University}
{    P.O. Box 7907}
{    Raleigh, NC 27695}

interface


	const
		TNColumnSize = 30;
		TNRowSize = 30;

	type
		TNColumnVector = array[1..TNColumnSize] of EXTENDED;
		TNRowVector = array[1..TNRowSize] of EXTENDED;
		TNmatrix = array[1..TNColumnSize] of TNRowVector;
		TNSquareMatrix = array[1..TNRowSize] of TNRowVector;
		TNString40 = string[40];
		FitType = (expo, fourier, log, poly, power, user);


	procedure LeastSquares (nstandards: INTEGER; xdata, ydata: TNColumnVector;
										{Array[1..MaxV] of EXTENDED}
									ncoefficients: INTEGER; var solution: TNRowVector; var yfit, residuals: TNColumnVector;
										{Array[1..MaxV] of EXTENDED}
									var fitsd, variance: EXTENDED; var err: BYTE; typeoffit: FitType); {Poly,Expo,Power,Log}


implementation

	procedure LeastSquares;

		const
			MaxN = 30;

		var
			i, j, k, L, nc, nx, ct: INTEGER;
			a: array[1..6, 1..7] of EXTENDED;
			s, t, valu: array[1..7] of EXTENDED;
			z, p, r, tval, temp: EXTENDED;
			yv, xv: array[1..MaxN] of EXTENDED;
			xa: array[1..5, 1..MaxN] of EXTENDED;

	begin
		{check limits}
		nc := ncoefficients - 1;

		if (nstandards < 2) then begin
				err := 1;
				EXIT(LeastSquares);
			end;
		if (nc < 1) or (nc > 5) then begin
				err := 2;
				EXIT(LeastSquares);
			end;
		if (nc >= nstandards) then begin
				err := 3;
				EXIT(LeastSquares);
			end;

		ct := nstandards;
		if (ct > MaxN) then
			ct := MaxN;

		{initialize arrays}
		for j := 1 to 6 do
			for k := 1 to 7 do
				a[j, k] := 0;

		for j := 1 to 7 do begin
				s[j] := 0;
				t[j] := 0;
				valu[j] := 0;
			end;

		{copy incoming data and prepare columns of values}
		for i := 1 to ct do begin
				xv[i] := xdata[i];
				yv[i] := ydata[i];
				solution[i] := 0;
				yfit[i] := 0;
				residuals[i] := 0;
			end;
		fitsd := 0;

		{if type of fit is poly, create x^2, x^3, ...}
		if (typeoffit = poly) then
			if nc > 1 then
				for i := 1 to ct do begin
						tval := xv[i];
						temp := tval;
						for j := 2 to nc do begin
								tval := tval * temp;
								xa[j, i] := tval;
							end;
					end;

		{if type of fit is log, y=a ln(bx) then create ln(x)}
{		 - solve as y= a*ln(b) + a*ln(x)}
		if typeoffit = log then
			for i := 1 to ct do begin
					tval := xv[i];
					if tval > 0 then
						xv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
				end;

		{if type of fit is power, y=ax^b then create ln(y) and ln(x)}
{		  - solve as ln(y)=ln(a)+b*ln(x)}
		if typeoffit = power then
			for i := 1 to ct do begin
					tval := xv[i];
					if tval > 0 then
						xv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
					tval := yv[i];
					if tval > 0 then
						yv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
				end;

		{if type of fit is expo, y=ae^bx then create ln(y)}
{		 - solve as ln(y)=ln(a)+bx}
		if typeoffit = expo then
			for i := 1 to ct do begin
					tval := yv[i];
					if tval > 0 then
						yv[i] := LN(tval)
					else begin
							err := 4;
							EXIT(LeastSquares);
						end;
				end;

		{Build the matrix}
		nx := ncoefficients;
		for i := 1 to ct do begin
				valu[1] := 1.0;
				valu[2] := xv[i];
				if (typeoffit = poly) and (nc > 1) then
					for j := 2 to nc do
						valu[j + 1] := xa[j, i];
							{subscript order?}

				valu[nx + 1] := yv[i];
				for k := 1 to nx do
					for L := 1 to nx + 1 do begin
							a[k, L] := a[k, L] + valu[k] * valu[L];
							s[k] := a[k, nx + 1];
						end;

				s[nx + 1] := s[nx + 1] + valu[nx + 1] * valu[nx + 1];
			end; {for i}

		{Solve the matrix by inversion}
		for i := 2 to nx do
			t[i] := a[1, i];
		for i := 1 to nx do begin
				j := i - 1;
				repeat
					j := j + 1;
				until (j > nx) or (a[j, i] <> 0);

				if j > nx then begin
						err := 4;
						EXIT(LeastSquares);
					end;

				for k := 1 to nx + 1 do begin
						p := a[i, k];
						a[i, k] := a[j, k];
						a[j, k] := p;
					end;

				z := 1 / a[i, i];
				for k := 1 to nx + 1 do
					a[i, k] := z * a[i, k];

				for j := 1 to nx do
					if j <> i then begin
							z := -a[j, i];
							for k := 1 to nx + 1 do
								a[j, k] := a[j, k] + z * a[i, k];
						end;
			end; {for i=1..nx}

		{copy back the solution vector}
		for j := 1 to nx do
			solution[j] := a[j, nx + 1];

		if typeoffit = log then begin
				if solution[2] <= 0 then begin
						err := 4;
						EXIT(LeastSquares);
					end
				else begin
						p := solution[2];
						solution[2] := EXP(solution[1] / p);
						solution[1] := p;
					end;
			end;
		{solved as y= a*ln(b)+a*ln(x)}

		if typeoffit = power then
			solution[1] := EXP(solution[1]);
		{solved as ln(y)=ln(a)+b*ln(x)}

		if typeoffit = expo then
			solution[1] := EXP(solution[1]);
		{solved as ln(y)=ln(a)+b*ln(x)}

		{Create the predicted and residual columns}
		for i := 1 to ct do begin
				valu[1] := xv[i];
				if nc > 1 then
					for j := 2 to nc do
						valu[j] := xa[j, i];

				valu[nx] := yv[i];
				p := solution[1];
				for j := 1 to nx - 1 do
					p := p + solution[j + 1] * valu[j];

				if (typeoffit = power) or (typeoffit = expo) then
					p := EXP(p); {convert from ln(y)}
				yfit[i] := p;
				residuals[i] := valu[nx] - p;
			end;
		p := 0;
		for i := 2 to nx do
			p := p + a[i, nx + 1] * (s[i] - t[i] * s[1] / ct);

		r := s[nx + 1] - s[1] * s[1] / ct;
		z := r - p;
		L := ct - nx;

		{Coeff. of determination (R-squared)= p/r}
{		 Std error of estimate=Sqrt(Abs(Z/L))}
		variance := z / L;
		fitsd := SQRT(ABS(z / L));
		err := 0;
	end; {LeastSquares}

end. {unit leastsq}