function [W,pvalue]=ranktest(Y,T);

% calculates Wilcoxon statistic and simulates one-sided p-value 

n=length(Y);
R=zeros(n,1);

sorty=sort(Y);

r=1;
while r<(n+1),
  yindex_r = find(Y==sorty(r));
  nr = length(yindex_r);
  R(yindex_r) = mean(r:r+nr-1);
  r = r+nr;
end;

% wilcoxon rank statistic
W = sum(R(T==1));

% simulation of p-value

J=100000;
Wdraws = zeros(J,1);
for j=1:J,
   Z = rand(n,1);
   [zs, Rdraw] = sort(Z);
   Wdraws(j) = sum(Rdraw(T==1));
end;    

% one-sided p-value
pvalue = sum(Wdraws>W)/J;