procedure out(t:matrix);
var
i1,j1:integer;
begin
for i1:=1 to m do
begin
for j1:=1 to m do
write(' ',t[i1,j1]:3:3,' ');
writeln
end
end;
{===================================================================}
Begin
writeln('ПРОГРАММА РАСЧЕТА ГЛАВНЫХ КОМПОНЕНТ ПО ЗАДАННОМУ РАСПРЕДЕЛЕНИЮ');
writeln;
x[1,1]:=2;x[1,2]:=1.3;x[1,3]:=0.55;x[2,1]:=4;x[2,2]:=1.42;x[2,3]:=5.1
x[3,1]:=1.1;x[3,2]:=5.3;x[3,3]:=0.55;x[4,1]:=2.14;x[4,2]:=5.12;x[4,3]:=1.9;
{------стандартизуем значения признаков-----------}
for j:=1 to m do
begin
{----находим среднее и сигму-----}
s:=0;x_:=0;
for i:=1 to n do
s:=s+x[i,j];
x_:=s/n;s:=0;
for i:=1 to n do
s:=s+(x[i,j]-x_)*(x[i,j]-x_);
s:=sqrt(s/n);
{------нормируем-------}
for i:=1 to n do
z[i,j]:=(x[i,j]-x_)/s
end;
{---------находим матрицу парных корреляций R=(1/n)*Z'*Z----------}
for j:=1 to m do
for i:=1 to m do
begin
s:=0;
for k:=1 to n do
s:=s+z[k,j]*z[k,i];
r[j,i]:=s/n
end;
{-------------выводим матрицу R------------}
writeln('Матрица парных корреляций R:');
out(r);
{-------=====находим собственные числа матрицы R======----------}
{-----приравниваем R и _a_-------}
for i:=1 to m do
for j:= 1 to m do
_a_[i,j]:=r[i,j];
p[1]:=3;{т.к на главной диагонали единицы}
for i:=1 to m do
for j:=1 to m do
if i<>j
then
_b_[i,j]:=_a_[i,j]
else
_b_[i,j]:=-2;
for q:=2 to m do
{----вычисляем p[q] и определитель-----}
begin
{----вычисляем A[q]----}
for i:=1 to m do
for j:=1 to m do
begin
s:=0;
for k:= 1 to m do
s:=s+r[i,k]*_b_[k,j];
a_[i,j]:=s
end;
{------вычисляем p[q]-------}
s:=0;
for i:=1 to m do
s:=s+a_[i,i];
p[q]:=s/q;
{----вычисляем B[q]-----}
for i:=1 to m do
for j:=1 to m do
if i<>j
then
b_[i,j]:=a_[i,j]
else
b_[i,j]:=a_[i,j]-p[q];
{----присваиваем предыдущим переменным значения текущих-----}
for i:= 1 to m do
for j:=1 to m do
begin
_a_[i,j]:=a_[i,j];
_b_[i,j]:=b_[i,j]
end
end;
{---------===решаем характеристическое уравнение===----------}
p[0]:=1;
for i:=1 to m do
p[i]:=-p[i];
for i:=1 to m do
for j:=1 to m do
l[i,j]:=0;
{------задаем начальные приближения------}
for i:=1 to m do
l[i,i]:=-p[i]/p[i-1];
{------выполняем итерационный процесс по методу Ньютона--------}
repeat
w:=0;
for i:=1 to m do
begin
b:=0;_b:=0;
{-----вычисляем значение полинома в i-й точке-------}
for j:=0 to m do
begin
s:=1;
for k:=0 to m-j-1 do
s:=s*l[i,i];
b:=b+p[j]*s
end;
{------находим максимальную невязку-------}
if b>w then
w:=b;
{------вычисляем значение производной в i-й точке------}
for j:=0 to m-1 do
begin
s:=1;
for k:=0 to m-j-2 do
s:=s*l[i,i];
_b:=_b+(m-j)*p[j]*s
end;
{------вносим поправку для i-й точки-------}
l[i,i]:=l[i,i]-(b/_b)
end
{----выходим из процесса при достижении требуемой точности----}
until w<0.0001;
{-------выводим собственные числа на экран---------}
writeln('Собственные числа матрицы R:');
for i:=1 to m do
writeln('L[',i,'] := ',l[i,i]:3:3);
{-----======находим матрицу собственных векторов u======---------}
{-----последним компонентам придаем единичные значения-----}
for i:= 1 to m do
u[m,i]:=1;
{------==решаем m систем уравнений==------}
for q:=1 to m do
begin
{----заполняем левые части-----}
for i:=1 to m-1 do
for j:=1 to m-1 do
if i=j
then
c[i,j]:=1-l[q,q]
else
c[i,j]:=r[i,j];
{----заполняем правые части-----}
for i:=1 to m-1 do
d[i]:=-r[i,m]*u[m,i];
{---------решаем систему методом Гаусса-----------}
i:=1;
{-------------прямой ход---------------}
repeat
{---нормируем элементы i-й строки---}
d[i]:=d[i]/c[i,i];
for j:=m-1 downto i do
c[i,j]:=c[i,j]/c[i,i];
{----делаем нули под ведущим элементом----}
for k:=i+1 to m-1 do
begin
d[k]:=d[k]-d[i]*c[k,i];
for j:=m-1 downto i do
c[k,j]:=c[k,j]-c[i,j]*c[k,i]
end;
i:=i+1
until i=m;
{------------обратный ход-------------}
u[m-1,q]:=d[m-1];
for i:=m-2 downto 1 do
begin
u[i,q]:=d[i];
for j:=i+1 to m-1 do
u[i,q]:=u[i,q]-u[j,q]*c[i,j]
end
end;
{------нормируем собственные векторы - находим матрицу v---------}
for j:=1 to m do
begin
s:=0;
for i:=1 to m do
s:=s+u[i,j]*u[i,j];
for i:=1 to m do
v[i,j]:=u[i,j]/sqrt(s)
end;
{--выводим нормированную матрицу собственных векторов на экран---}
writeln('Матрица нормированных собственных векторов V:');
out(v);
{---------находим матрицу факторного отображения a----------}
for i:=1 to m do
for j:=1 to m do
begin
s:=0;
for k:=1 to m do
s:=s+v[i,k]*sqrt(l[k,j]);
a[i,j]:=s
end;
{--------выводим матрицу факторного отображения---------}
writeln('Матрица факторного отображения A:');
out(a);
{===находим матрицу, обратную a, методом m-кратного пересчета===}
for k:=1 to m do
{-----цикл пересчета-----}
begin
for i:=1 to m do
for j:=1 to m do
if (i=k) or (j=k)
then
if i=j
then
a_1[i,j]:=-1/a[i,j]
else
a_1[i,j]:=-a[i,j]/a[k,k]
else
a_1[i,j]:=a[i,j]-a[i,k]*a[k,j]/a[k,k];
for i:=1 to m do
for j:=1 to m do
a[i,j]:=a_1[i,j]
end;
for i:=1 to m do
for j:=1 to m do
a_1[i,j]:=-a[i,j];
{------===находим матрицу значений главных компонент F===-------}
for i:=1 to m do
for j:=1 to n do
begin
s:=0;
for k:=1 to m do
s:=s+a_1[i,k]*z[j,k];
f[i,j]:=s
end;
{-------выводим матрицу F на экран в транспонированном виде-------}
writeln('Матрица значений главных компонент F:');
for i:=1 to n do
begin
for j:=1 to m do
write(' ',f[j,i]:3:3,' ');
writeln
end;
writeln;
readln;
End.