function g_perm(z,W,n)

%G_PERM performs a random-permutation test of G and G* 
%stats for a given set of area-unit data. Here G* is defined 
%to be
%
%     G* = (z'*W*z)/(z'*U*U'*z)
%
%for ANY nonnegative weight matrix, W,(usually with nonzero 
%diagonal elements), where U is a unit column vector. G is then 
%defined to be the same expression with diagonal element 
%removed from W.

%Written by: TONY E. SMITH, 4/1/98


%INPUTS: (i)    z = Nx1 vector of data values for test
%        (ii)   W = NxN weight matrix of spatial relations 
%        (iii)  n = number of simulations to be run
%
%SCREEN OUTPUT: ROW 1 = value and significance level of G.
%                       (example: OUT = .1 implies that there 
%                       is only a 10% chance of realizing a 
%                       G value as high as the observed 
%                       value in a random permutation of the
%                       given values over locations.)
%               ROW 2 = same for G*.


N = length(z);

U = ones(N,1) ;

M = zeros(n,2) ;

OUT = zeros(2,2) ;

%Compute G and G*(= GG)


GG = (z'*W*z)/(z'*U*U'*z) ;

W0 = W - diag( diag( W ) ) ; %This removes diagonals

G = (z'*W0*z)/(z'*U*U'*z) ;

OUT(1,1) = G ;
OUT(2,1) = GG ;


%Compute test values

for i = 1:n
   
   %Compute random permutation of locations
   
   P = randperm(N);
   
  %Compute G and G* (= GG)


	GG = (z(P)'*W*z(P))/(z(P)'*U*U'*z(P)) ;

	G = (z(P)'*W0*z(P))/(z(P)'*U*U'*z(P)) ;
   
      
   M(i,:) = [G,GG] ;
   
end

%Rank test values for I,rho,r (assumes nonnegative correlation)

for j=1:2 
   
   v = - sort(-M(:,j)); %Sorts in descending order

	del = 1 ;

	i = 1 ;
   
   prob = 1 ;
	
	while del > 0
   
   	if OUT(j,1) >= v(i)
      
      	prob = i/(n+1) ;
      
      	del = 0;
      
   	end
   
   	i = i + 1 ;
   
   end
   
   OUT(j,2) = prob ;
   
end

%Prepare Screen Output

TABLE(1,:) = {'INDEX','VALUE','PROB'}; %Column labels

TABLE(2:3,1) = {'G','G*'};  %Row labels

TABLE(2:3,2:3) = num2cell(OUT);

CLUSTERING_RESULTS = TABLE

