%
%  Calculates the inverse normal distribution.  
%
%   x=phiinv(z) 
%
%  gives z=int((1/sqrt(2*pi))*exp(-t^2/2),t=-infinity..x)
%
function x=phiinv(z)
if (z >= 0.5), 
  x=sqrt(2)*erfinv((z-0.5)/.5);
else
  x=-phiinv(1-z);
end