Inverse of the Standard Normal Cumulative Distribution

Anyone have a function akin to Excel's NORM.S.INV?

It returns the inverse of the standard normal cumulative distribution. That is, it returns the Z value, given a supplied probability value. The standard distribution has a mean of zero and a standard deviation of one. Excel uses iteration over NORM.S.DIST to find the Z value that matches.

You could save us some time and share the complete formula if its available otherwise you'd have to wait for someone who replicated it

I don't really care about a specific Z value or probability, so I'm wondering if this is equivalent:

1 / NORMDIST(mtrandom(),0,1)

I'm trying to replicate this formula:

random value = standard deviation * NORMSINV(RAND())

There is no closed form for the inverse of the CDF for a normal distribution.

Your formula is not correct (1 / NORMDIST(mtrandom(),0,1))

There are good approximations that you can use - Please refer to this link:

Inverse Normal CDF approximations

FWIW, I've once made two AFL versions of NormSInv:

1st version applies Abramowitz and Stegun approximation. I've named it asNormSInv ('as...' stands for Abramowitz and Stegun). (A Python version may be found here.)

function asNormSInv( p ) {
	/// posted at official AmiBroker forum
	/// @link https://forum.amibroker.com/t/inverse-of-the-standard-normal-cumulative-distribution/10072/5
	/// applies Abramowitz and Stegun approximation  
	/// at 26.2.23. of "Handbook of Mathematical Functions"
	/// AFL version by fxshrat@gmail.com, AB 6.10+ required
	local j, t, m, ra, result;
	if ( typeof(p) != "number" )
		Error("asNormSInv(p): 'p' must be type 'number'!");	
	if (p <= 0 || p >= 1) 
        Error(StrFormat("asNormSInv(p): Invalid argument! p(%g) must be larger than 0 but less than 1.", p)); 
	/// 
	m = MxFromString( "[[2.515517,1.432788],[8.02853e-1,1.89269e-1],[1.0328e-2,1.308e-3],[0,0]]" ); 
	/// 
	t = sqrt(-2*log(IIf(p < 0.5, p, 1-p)));    
	for( j = 0; j < 2; j++ )	m[3][j] = (m[2][j]*t + m[1][j])*t + m[0][j];// y = a*x^2+b*x+c    
	ra = t - m[3][0] / (m[3][1]*t+1); 
	///
	return IIf(p < 0.5, -ra, ra);
}

Sample output

EnableTextOutput(0);
str = "[0.000001;
        0.00001;
        0.001;
        0.05;
        0.15;
        0.25;
        0.35;
        0.45;
        0.55;
        0.65;
        0.75;
        0.85;
        0.95;
        0.999;
        0.99999;
        0.999999]";

mat = MxFromString(str);
printf(MxToString(mat));

printf( "\n\nNormSInv(p):" );
rownum = MxGetSize( mat, 0 );
for( i = 0; i < rownum; i++ ) {
	p = mat[i][0];
	printf( "\n%1.4f", asNormSInv(p) );
}

AB output of 1st version

24


2nd AFL version applies Peter J. Acklam approximation. I've named it paNormSInv ('pa...' stands for Peter Acklam). As for polynomial evaluation I have made Horner function (-> Horner's rule).

function MxHornerValue(mat, x) {		
	/// applies Horner's rule for polynomial evaluation
	/// https://en.wikipedia.org/wiki/Horner%27s_method
	/// AFL version by fxshrat@gmail.com, AB 6.00+ required	
	local cs, i, rownum;
	if ( typeof(mat) != "matrix" )
		Error("MxHornerValue(): 1st argument must be type 'matrix'!");
	if ( typeof(x) != "number" )
		Error("MxHornerValue(): 2nd argument must be type 'number'!" );	
	cs = 0; 
	rownum = MxGetSize(mat, 0);
	for (i = rownum-1; i >= 0; i--)
		cs = cs * x + mat[i][0];
	return cs;
}

function paNormSInv(p) {
	/// posted at official AmiBroker forum
	/// @link https://forum.amibroker.com/t/inverse-of-the-standard-normal-cumulative-distribution/10072/5
	/// function applies Peter J. Acklam approximation to
	/// compute the inverse normal cumulative distribution. 
	/// AFL version by fxshrat@gmail.com, AB 6.10+ required	
	local p_low, p_high, q, t, a, b, ra;
	if ( typeof(p) != "number" )
		Error("paNormSInv(p): 'p' must be type 'number'!" );	
	if ( p <= 0 || p >= 1 ) 
		Error(StrFormat("paNormSInv(p): Invalid argument! p(%g) must be larger than 0 but less than 1.", p)); 
	//
	p_low = 0.02425; 
	p_high = 1 - p_low;
	//
	if ( p >= p_low && p <= p_high ) {
		q = p - 0.5; 
		t = q * q;
		a = MxFromString("[2.506628277459239;-30.66479806614716;138.3577518672690;-275.9285104469687;220.9460984245205;-39.69683028665376]");	
		b = MxFromString("[-13.28068155288572;66.80131188771972;-155.6989798598866;161.5858368580409;-54.47609879822406]");	
	} else {	
		q = 1;
		t = sqrt(-2*log(IIf(p > p_high, 1-p, p)));	
		a = MxFromString("[2.938163982698783;4.374664141464968;-2.549732539343734;-2.400758277161838;-3.223964580411365e-1;-7.784894002430293e-3]");		
		b = MxFromString("[3.754408661907416;2.445134137142996;3.224671290700398e-1;7.784695709041462e-3]");		
	}
	//
	ra = MxHornerValue(a, t) * q / (MxHornerValue(b, t) * t + 1);	
	return IIf(p > p_high, -ra, ra);
}

AB sample output 2nd version:

EnableTextOutput(0);
str = "[0.000001;
        0.00001;
        0.001;
        0.05;
        0.15;
        0.25;
        0.35;
        0.45;
        0.55;
        0.65;
        0.75;
        0.85;
        0.95;
        0.999;
        0.99999;
        0.999999]";

mat = MxFromString(str);
printf(MxToString(mat));

printf( "\n\nNormSInv(p):" );
rownum = MxGetSize( mat, 0 );
for( i = 0; i < rownum; i++ ) {
	p = mat[i][0];
	printf( "\n%1.4f", paNormSInv(p) );
}

39


Excel output for comparison:

22

11 Likes

You can use the same as Excel uses - iterative approach. Since CDF is monotonically increasing a very fast binary search can be used with built-in NormDist function

// iterative approach (binary search)

function InvNormDist( y )
{
  if( y > 0 AND y < 1 )
  {
	  x = 0;
	  step = 6;
	  
	  do
	  {
		y1 = NormDist( x ); // the function to calculate inverse of

		if( y1 > y ) 
			x -= step;
		else 
			x += step;
		
		step *= 0.5;
	  }
	  while( step > 1e-9 );
  }
  else
  {
    //Error("Out of range");
	x = Null;
  }
  
  
  return x;
}

for( i = 1; i < 100; i++ )
{
  y = i/100;
  printf("%g -> %g\n", y, InvNormDist( y ) );
}

The beauty of binary search is that it can be used to quickly find the inverse of any monotonic function by replacing single line of code.

15 Likes