function g2_rand(N1,N2,W,sims)

%G2_RAND tests G* for pop1 versus a bigger pop2 by subsampling pops of 
%size pop1 from pop2, and comparing G*(1) with these samples. 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. 

%Written by: TONY E. SMITH, 4/30/03


%INPUTS: (i)    N1   = nx1 vector of pop1 values for n regions
%        (ii)   N2   = nx1 vector of pop2 values for n regions
%        (iii)  W    = weight matrix of spatial relations 
%        (iii)  sims = number of simulations to be run
%
%SCREEN OUTPUT: P-value for G*(1) relative to sample values G*(i), i=1,..,sims

n1 = length(N1);

U = ones(n1,1) ;

%Compute G* = G


G = (N1'*W*N1)/(N1'*U*U'*N1) ;


%Construct Sampling Probabilities

Prob = N2/sum(N2);

cdf = cumsum(Prob);

GG = zeros(sims,1);

for s = 1:sims
   
   %Construct Sample G* values
   
   s1 = sum(N1);
   
   S = zeros(n1,1);
   
   for i = 1:s1
      
      r = rand(1,1);
      
      index = 1;
      
      j = 1;
      
      while index > 0
         
         if r <= cdf(j)
            
            S(j) = S(j) + 1;
            
            index = 0;
            
         else
            
            j = j + 1;
            
         end
         
      end %End while
      
   end %End For
   
   GG(s) = (S'*W*S)/(S'*U*U'*S);
   
end %End construction


%Rank test values for GG


v = - sort(-GG); %Sorts in descending order

del = 1 ;

i = 1 ;

prob = 1 ;

while del > 0
   
   if i == sims
      
      prob = 1;
      
      del = 0;
      
   elseif G >= v(i)
            
      prob = i/(sims+1) ;
      
      del = 0;
      
   end
   
   i = i + 1 ;
   
end


%Prepare Screen Output

disp(['P-Value = ',num2str(prob)]);