Смекни!
smekni.com

Математическая модель метода главных компонент (стр. 2 из 2)

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.