Writing a summation code where i≠j

2 views (last 30 days)
aposta95
aposta95 on 10 Jun 2022
Commented: Jan on 10 Jun 2022
Hi, I want to write a summation code to calculate Getis ord index as below:
.
Here, x is a 1xn vector, and w is a square matrix where w_ii=0
I wrote the code to calculate the denominator, but it's quite complicated to calculate the nominator.
I tried with the below code, but I'm not sure it's right...
How can I do that?
x = [1,2,3,4]; % feature vector
w= rand(4,2); D=pdist(w);
W= squareform(D); % weight matrix
% Calculate Denominator
denom=0;
for i = 1:4
for j = 1:4
if i~=j
denom = denom + x(j)*x(i);
end
end
end
% Calculate nominator
Wu = triu(W,1);
Wd = tril(W,-1);
nom = (sum(sum(Wu)) + sum(sum(Wd))) * denom;

Accepted Answer

Voss
Voss on 10 Jun 2022
The weight matrix is zeros along the diagonal. Thus you can use it as-is in the summation, since the terms where i==j (which correspond to the diagonal elements of W) will be 0.
To get the same effect in the denominator, you can set the diagonal elements of x.*x.' to zero and use that.
x = [1,2,3,4]; % feature vector
w = rand(4,2);
D = pdist(w);
W = squareform(D) % weight matrix
W = 4×4
0 0.4616 0.3697 0.5919 0.4616 0 0.8026 0.7224 0.3697 0.8026 0 0.8581 0.5919 0.7224 0.8581 0
denom_matrix = x.*x.';
denom_matrix(1:numel(x)+1:end) = 0
denom_matrix = 4×4
0 2 3 4 2 0 6 8 3 6 0 12 4 8 12 0
G = sum(W.*x.*x.','all')/sum(denom_matrix,'all')
G = 0.7226
For reference, here is a way to do the same thing using for loops:
denominator = 0;
numerator = 0;
for ii = 1:numel(x)
for jj = 1:numel(x)
if ii == jj % skip the ii == jj terms
continue
end
denominator = denominator + x(ii)*x(jj);
numerator = numerator + W(ii,jj)*x(ii)*x(jj);
end
end
G_loop = numerator/denominator
G_loop = 0.7226
Check that the results are identical:
isequal(G,G_loop)
ans = logical
1
  2 Comments
aposta95
aposta95 on 10 Jun 2022
Many thanks for your elaboration :)

Sign in to comment.

More Answers (2)

Abhishek Chakram
Abhishek Chakram on 10 Jun 2022
Hi,
As I understood, you want to calculate the numerator of Getis ord index.
You can do it the same way as you calculated the denominator. The code will look like :
numerator = 0;
for i = 1:4
for j = 1:4
if i~=j
numerator = numerator + W(i,j)*x(i)*x(j);
end
end
end

David Hill
David Hill on 10 Jun 2022
No loops needed
G=sum(w.*x.*x','all')/sum(x.*x','all');

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!