icdf versus erfinv?

What is the difference between using the icdf function and the erfinv function to compute the inverse CDF for a normal distribution? Consisder the example below. When I choose a value for p in the 'heart' of the distribution (e.g., p = 0.2) or in the upper tail (e.g., p = 1 - 1e-9), then the difference between x1 and x2 is zero. However when I choose a value for p close to zero, as shown in the example, there is a non-zero difference between x1 and x2. Which of these two computations is more accurate, and why?

% Mean and standard deviation of a normal distribution
mu  = 3.0;
sig = 1.5;

% Quantile for which we want to find the associated value of x
p = 1e-9;

% Use the icdf function to determing x at p
x1 = icdf('Normal',p,mu,sig);

% Now define an anonymous function to perform the same computation, using the erfinv function
myinversefun = @(p,m,s) (s*sqrt(2))*erfinv(2*p-1) + m;
x2 = myinversefun(p,mu,sig);

% Print the difference in results in exponential format
fprintf('difference = %16.9e\n',x1-x2);

Answer 1

[John D'Errico]

In general, use the ICDF code!

Why would it be more accurate? Because it can be more intelligent, using a more accurate approximation in the tails. For example, using erfinv as you did does something that is dangerous in the tails.

x1
x1 =
         -5.99671052251153

x2
x2 =
          -5.9967105158771

vpa((sig*sqrt(2))*erfinv(2*sym(p)-1) + mu)
ans =
-5.9967105225115303073434653073673

So x2 is your computation. Then I computed a symbolic version of your computation. As you can see, it is now the same as x1.

The point is, icdf gave you a result that is correct to the last digit printed, whereas the double precision erfinv code failed.

[John Schaefer]

Thanks for the explanation, John! The example with vpa is instructive.