function x = norminv(p)
% x = norminv(p) calculates the inverse of the normal 
%        distribution function, i.e. the x such that 
%        P(X <= x) = p, where X is N(0,1) random variable.
%      
% Uses Newton's method employing NORMDENS and NORMDIST.
%

% Written by Ulf Brannlund, November 9, 1997

if min(size(p)) ~= 1; error('Input p needs to be scalar.'); end	  
if ~isreal(p) ; error('Input p needs to be real.'); end	  
if p < eps ; error('Input p needs to be positive'); end	  
if p > 1-eps ; error('Input p needs to be less than 1'); end	  

x0 = 0;
x1 = inf;
while abs(x0-x1) > 1e-14
  x1 = x0;
  x0 = x0-(normdist(x0)-p)/normdens(x0);
end
x = x0;
    
