function x = norminv(p,mu,sigma);
%NORMINV Inverse of the normal cumulative distribution function (cdf).
%	X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
%	mean, MU, and standard deviation, SIGMA.
%
%	The size of X is the common size of the input arguments. A scalar input  
%	functions as a constant matrix of the same size as the other inputs.	 
%
%	Default values for MU and SIGMA are 0 and 1 respectively.

%	References:
%	   [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%	   Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2

%	Copyright (c) 1993 by The MathWorks, Inc.
%	$Revision: 1.1 $  $Date: 1993/05/24 18:55:33 $

if nargin < 3, 
    sigma = 1;
end

if nargin < 2;
    mu = 0;
end

[errorcode p mu sigma] = distchck(3,p,mu,sigma);

if errorcode > 0
    error('The arguments must be the same size or be scalars.');
end

% Allocate space for x.
x = zeros(size(p));

% Return NaN if the arguments are outside their respective limits.
k = find(sigma <= 0 | p < 0 | p > 1);
if any(k)
    x(k) = NaN * ones(size(k)); 
end

% Put in the correct values when P is either 0 or 1.
k = find(p == 0);
if any(k)
    x(k) = -Inf * ones(size(k));
end

k = find(p == 1);
if any(k)
    x(k) = Inf * ones(size(k)); 
end

% Compute the inverse function for the intermediate values.
k = find(p > 0  &  p < 1 & sigma > 0);
if any(k),
    x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k);
end
