I'm not familiar with the algoritm, so it could still have some mistakes left...
However, i was able to identify a couple of mistakes:
- x1 should be indexed with i1 and x0 with i0
- you first need to update j before you can use it as an index into x2
- the condition states that you should enter if i2 is not equal to i1
- i'm guessing that the matrix A you provided does not have the correct dimensions, i assume it should be a 7x7 matrix?
- also min( median( ) ) is not the same as the low median, since x2 is a vector the median will return a scalar
Below i used an random A matrix, hence the 'PS' result is not the same as the one you expect. Where did you obtain the matrix A and the expected result?
A = [10 -10; 1 0; -1 0; 0 -1 ; 0 1; 11 -10; -1 -1];
% use a random A matrix...
A = rand(7) * 100;
C = cov(A)^(-1/2) * A;
M = size(A,1);
omega = C*C';
x0 = zeros(M ,1);
x1 = zeros(M ,1);
x2 = zeros(M-1,1); % 1 element less since we neglect the case where i2 == i1
PS = zeros(M ,1);
for i0 = 1:M
for i1 = 1:M
j = 0;
for i2 = 1:M
if i2 ~= i1
j = j + 1;
x2( j ) = abs( omega(i1,i0) + omega(i2,i0) );
end
end
x1( i1 ) = median(x2,'all');
end
x0( i0 ) = 1.1 * median(x1,'all');
end
for i3 = 1:M
PS(i3) = max( omega(i3)/x0 );
end
PS
