icdf versus erfinv?

5 vues (au cours des 30 derniers jours)
John Schaefer
John Schaefer le 3 Oct 2019
Commenté : John Schaefer le 7 Oct 2019
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);

Réponse acceptée

John D'Errico
John D'Errico le 3 Oct 2019
Modifié(e) : John D'Errico le 3 Oct 2019
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.
  1 commentaire
John Schaefer
John Schaefer le 7 Oct 2019
Thanks for the explanation, John! The example with vpa is instructive.

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by