cdfnor

// Compute p from x
x = -1;
Mean = 0;
Std = 1;
// P = 0.1586553, Q = 0.8413447
[P,Q]=cdfnor("PQ",x,Mean,Std)

// Compute p from a very rare x
// Notice that Q should be used when P is close to 1, 
// while P should be used when Q is close to 1.
x = 10;
Mean = 0;
Std = 1;
// P = 1, Q = 7.620D-24
[P,Q]=cdfnor("PQ",x,Mean,Std)

// Make a table of P,Q
Mean = 0;
Std = 1;
x = -6:6;
[P,Q]=cdfnor("PQ",x,Mean*ones(x),Std*ones(x));
[x' P' Q']

// Find x from p,q
Mean = 0;
Std = 1;
p = 0.1;
q = 0.9;
// x = - 1.2815516 
x = cdfnor("X",Mean,Std,p,q)

// Find Mean from p,q
x = 2;
Std = 1;
p = 0.1;
q = 0.9;
// Mean = 3.2815516
Mean = cdfnor("Mean",Std,p,q,x)

// Find Standard Deviation from p,q
Mean = 0;
p = 0.1;
q = 0.9;
x = 2;
// Std = - 1.5606083
Std = cdfnor("Std",p,q,x,Mean)

// Plot the function
h = scf();
Mean = [0 0 0 -2];
Std2 = [0.2 1.0 5.0 0.5];
cols = [1 2 3 4];
nf = size(cols,"*");
lgd = [];
for k = 1 : nf
  x = linspace(-5,5,1000);
  P=cdfnor("PQ",x,Mean(k)*ones(x),sqrt(Std2(k))*ones(x));
  plot(x,P)
  str = msprintf("Mean=%s, Std^2=%s",string(Mean(k)),string(Std2(k)));
  lgd($+1) = str;
end
for k = 1 : nf
h.children.children.children(nf - k + 1).foreground = cols(k);
end
legend(lgd);

// Plot the Inverse Cumulated Standard function
Mean = 0;
Std = 1;
p = linspace(1e-10,1-1e-10,1000);
q = 1-p;
x = cdfnor("X",Mean*ones(p),Std*ones(p),p,q);
plot(p,x)
xtitle("Inverse Cumulated Distribution Normal Standard Function","p","x");

// The Inverse Normal CDF is ill-conditionned when 
// p is close to p=0.5. 
// This is because, at p=0.5, the first derivative is non-zero,
// while the function is zero.
// This may produce a large relative error in the function value, 
// that is, some digits may not be accurate when p is getting 
// closer and closer to p=0.5.
Mean = 0;
Std = 1;
p = 0.500000001;
q = 1-p;
x = cdfnor("X",Mean,Std,p,q)
// The expected probability was computed from a symbolic software, 
// with the expression sqrt(2)*erfinv(2*p-1)
// where p = 1/2+10^-9.
expected = 2.50662827463100050e-9;
// The relative error is about 1.e-8, meaning that 
// only 8 digits are accurate.
abs(x-expected)/abs(expected)

// The zero-finder algorithm may be unable to find the standard 
// deviation in some situations when the inversion cannot be 
// performed.
Mean = 1;
p = 0.5;
q = 0.5;
x = 1;
// Std = Nan
Std = cdfnor("Std",p,q,x,Mean)