Метод наименьших квадратов

Очень востребованный численный метод обработки результатов экспериментальных измерений данных. Позволяет найти кривую, наилучшим образом описывающую полученные результаты. Программа составлена на языке Паскаль. Читайте комментарии в тексте.

program metod_naim_kwadr; 
{Метод наименьших квадратов} 
uses crt; 
const max=20; 
type massiv=array[1..3,1..4] of real; 
var 
xy:array[1..2,1..max] of real; 
b,a0,a1,a2,s:real; i,j,n:byte; 
procedure linapr(n:byte); 
var xsred,ysred,sum1,sum2,s:real; 
begin 
xsred:=0; ysred:=0; sum1:=0; sum2:=0; 
for i:=1 to n do 
begin 
xsred:=xsred+xy[1,i]; ysred:=ysred+xy[2,i]; 
sum1:=sum1+(xy[1,i]*xy[2,i]); 
sum2:=sum2+(xy[1,i]*xy[1,i]); 
end; 
xsred:=xsred/n;ysred:=ysred/n; 
sum1:=sum1/n; sum2:=sum2/n; 
b:=(sum1-(xsred*ysred)) / (sum2-(xsred*xsred)); 
a0:=ysred-(b*xsred); 
write('ПРЯМАЯ: '); writeln('y=',b:5:2,'* x + ',a0:5:2); 
s:=0; for i:=1 to n do s:=s+sqr(xy[2,i]-b*xy[1,i]-a0); 
writeln('Сумма квадратов отклонений = ',s:5:3); 
end; 
procedure stepenapr(n:byte); 
var xsred,ysred,sum1,sum2,s:real; 
lxy:array[1..2,1..max] of real; 
begin 
xsred:=0; ysred:=0; sum1:=0; sum2:=0; 
for i:=1 to 2 do for j:=1 to n do 
if xy[i,j]<>0 then lxy[i,j]:=ln(abs(xy[i,j])) 
else lxy[i,j]:=0; 
for i:=1 to n do 
begin 
xsred:=xsred+lxy[1,i]; ysred:=ysred+lxy[2,i]; 
sum1:=sum1+(lxy[1,i]*lxy[2,i]); 
sum2:=sum2+(lxy[1,i]*lxy[1,i]); 
end; 
xsred:=xsred/n;ysred:=ysred/n; 
sum1:=sum1/n; sum2:=sum2/n; 
b:=(sum1-(xsred*ysred)) / (sum2-(xsred*xsred)); 
a0:=ysred-(b*xsred); 
write('СТЕПЕННАЯ: '); writeln('y=',exp(a0):5:2,'* x ^ ',b:5:2); 
s:=0; for i:=1 to n do s:=s+sqr(xy[2,i]-exp(a0+b*lxy[1,i])); 
writeln('Сумма квадратов отклонений = ',s:5:3); 
end; 
procedure pokazatelnapr(n:byte); 
var xsred,ysred,sum1,sum2,s:real; 
ly:array[1..max] of real; 
begin 
xsred:=0; ysred:=0; sum1:=0; sum2:=0; 
for j:=1 to n do 
if xy[2,j]<>0 then ly[j]:=ln(abs(xy[2,j])) 
else ly[j]:=-1e10; 
for i:=1 to n do 
begin 
xsred:=xsred+xy[1,i]; ysred:=ysred+ly[i]; 
sum1:=sum1+(xy[1,i]*ly[i]); 
sum2:=sum2+(xy[1,i]*xy[1,i]); 
end; 
xsred:=xsred/n;ysred:=ysred/n; 
sum1:=sum1/n; sum2:=sum2/n; 
b:=(sum1-(xsred*ysred)) / (sum2-(xsred*xsred)); 
a0:=ysred-(b*xsred); 
write('ПОКАЗАТЕЛЬНАЯ: '); 
writeln('y=',exp(a0):5:2,'* ',exp(b):5:2,'^ x '); 
s:=0; for i:=1 to n do s:=s+sqr(xy[2,i]-exp(a0+b*xy[1,i])); 
writeln('Сумма квадратов отклонений = ',s:5:3); 
end; 
procedure giperbapr(n:byte); 
var xsred,ysred,sum1,sum2,s:real; 
z:array[1..max] of real; 
begin 
xsred:=0; ysred:=0; sum1:=0; sum2:=0; 
for j:=1 to n do 
if xy[1,j]<>0 then z[j]:=1/xy[1,j] 
else z[j]:=1e10; 
for i:=1 to n do 
begin 
xsred:=xsred+z[i]; ysred:=ysred+xy[2,i]; 
sum1:=sum1+(z[i]*xy[2,i]); sum2:=sum2+(z[i]*z[i]); 
end; 
xsred:=xsred/n;ysred:=ysred/n; sum1:=sum1/n; sum2:=sum2/n; 
b:=(sum1-(xsred*ysred)) / (sum2-(xsred*xsred)); 
a0:=ysred-(b*xsred); 
write('ГИПЕРБОЛА:'); writeln('y=',b:5:2,'/ x + ',a0:5:2); 
s:=0; for i:=1 to n do s:=s+sqr(xy[2,i]-(a0+b*z[i])); 
writeln('Сумма квадратов отклонений = ',s:5:3); 
end; 
function delta(a:massiv;j:byte):real; 
{Вычисление определителей системы уравненений; 
при j=0 - главный, при j=1,2,3 - по переменным} 
var i:byte; 
begin 
if j<>0 then for i:=1 to 3 do a[i,j]:=a[i,4]; 
delta:=a[1,1]*a[2,2]*a[3,3]+a[2,1]*a[3,2]*a[1,3]+ 
a[3,1]*a[1,2]*a[2,3]-a[1,1]*a[3,2]*a[2,3]- 
a[2,1]*a[1,2]*a[3,3]-a[3,1]*a[2,2]*a[1,3]; 
end; 
procedure kwadrapr(n:byte); 
var 
i,j:byte;a:massiv; delta0,x:real; 
begin 
for i:=1 to 3 do for j:=1 to 4 do a[i,j]:=0; 
for i:=1 to n do 
begin 
a[1,2]:=a[1,2]+xy[1,i];a[1,3]:=a[1,3]+sqr(xy[1,i]); 
a[1,4]:=a[1,4]+xy[2,i]; 
a[2,3]:=a[2,3]+sqr(xy[1,i])*xy[1,i] ; 
a[2,4]:=a[2,4]+xy[1,i]*xy[2,i] ; 
a[3,3]:=a[3,3]+sqr(sqr(xy[1,i])); 
a[3,4]:=a[3,4]+xy[2,i]*sqr(xy[1,i]); 
end; 
a[1,1]:=n; a[2,1]:=a[1,2]; a[2,2]:=a[1,3]; a[3,1]:=a[2,2]; 
a[3,2]:=a[2,3]; delta0:=delta(a,0); 
if delta0 <> 0 then 
begin 
a0:=delta(a,1)/delta0; a1:=delta(a,2)/delta0; 
a2:=delta(a,3)/delta0; 
end; 
write('ПАРАБОЛА: '); 
writeln('Y=',a2:6:2,'*x*x+ ',a1:6:2,'*x+ ',a0:6:2); 
s:=0; for i:=1 to n do 
begin 
x:=a0+a1*xy[1,i]+a2*sqr(xy[1,i]); s:=s+sqr(x-xy[2,i]); 
end; 
write('Cумма квадратов отклонений ');writeln('=',s:5:3); 
end; 
BEGIN 
clrscr; writeln('Уравнение регрессии (метод наименьших квадратов)'); 
repeat write('Число точек:'); read(n); until n>=3; 
writeln('Введите исходные данные по точкам, X и Y:'); 
for i:=1 to n do 
begin 
writeln('точка', i:2,':'); read(xy[1,i]); 
gotoXY(whereX+4,whereY-1);read(xy[2,i]); 
end; 
{n:=10; 
xy[1,1]:=0.01 xy[2,1]:=6.698 xy[1,2]:=0.02 xy[2,2]:=17.20 
xy[1,3]:=0.04 xy[2,3]:=28.64 xy[1,4]:=0.1 xy[2,4]:=45.45 
xy[1,5]:=0.2 xy[2,5]:=59.67 xy[1,6]:=0.3 xy[2,6]:=68.68 
xy[1,7]:=0.4 xy[2,7]:=75.42 xy[1,8]:=0.5 xy[2,8]:=80.86 
xy[1,9]:=0.6 xy[2,9]:=85.45 xy[1,10]:=0.7 xy[2,10]:=89.45} 
{xy[1,11]:= xy[2,11]:= xy[1,12]:= xy[2,12]:=}
linapr(n);kwadrapr(n);stepenapr(n);giperbapr(n);pokazatelnapr(n); 
repeat until keypressed; 
END.