% Point bi-serial rank-sum correlation function
%           __   __
%           Y1 - Y0
% r = 1 - -----------  SQRT (N1xN0/N(N-1))
%              Sy

% Each row of x should be a sample, and each column is a
% binary variable.  y is a non-binary variable represented
% as a column vector.

function r = biserial (x, y)

   % Get the sizes and build the result matrix

   [nr, nc] = size (x);
   

   r = zeros (nc);
   for v1 = 1:nc

      % Pull off the two vectors, and get the stats

      x1 = x(:, v1);
      x20 = y (find (x1==0));
      x21 = y (find (x1==1));
      M0 = mean (x20);
      M1 = mean (x21);
      C0 = max (size (x20));
      C1 = max (size (x21));
      sdev = std (y);
      nt = 2 * nr;
      r(v1) = (M1-M0)/sdev * sqrt ((C0*C1)/(nt*(nt-1)));
    
   end

