Java implementation of Excels statistical functions NORMINV

I was recently in need of a java implementation of the inverted normal distribution function, named “NORMINV” in Excel 2003 and later: NORMINV.

I only found a JavaScript implementation in kmpm’s repository norminv.js. It has quite a history, coming over C# from C++. I have the impression nobody wants to deal with the mathematics 😉

So here is my Java implementation of NORMINV. As an added bonus it implements org.apache.poi.ss.formula.functions.FreeRefFunction to be used as a User defined function in Apache POI:

import org.apache.poi.ss.formula.OperationEvaluationContext;
import org.apache.poi.ss.formula.eval.EvaluationException;
import org.apache.poi.ss.formula.eval.NumberEval;
import org.apache.poi.ss.formula.eval.OperandResolver;
import org.apache.poi.ss.formula.eval.ValueEval;
import org.apache.poi.ss.formula.functions.FreeRefFunction;
 
/**
 * Java Implementation of http://support.microsoft.com/kb/827358/en-us implementing 
 * Apache POIs {@link FreeRefFunction} to be used as a user defined function.
 * 
 * Mathematics found here: 
 * https://gist.github.com/kmpm/1211922/
 * https://gist.github.com/kmpm/1211922/raw/a11e0dfc9fab493bcdadc669f3213d11f1897ebf/norminv.js
 * 
 * Register for a given workbook with:
 * 
    <code>
		final UDFFinder udfs = new DefaultUDFFinder(new String[]{ "MS_NormInv" }, new FreeRefFunction[]{ new NormInv() }) ;        
		workbook.addToolPack(new AggregatingUDFFinder(new UDFFinder[] {udfs}));
    </code>
 * @author michael.simons
 */
public class NormInv implements FreeRefFunction {
	public ValueEval evaluate(ValueEval[] args, OperationEvaluationContext ec) {		
		try {
			final ValueEval p = OperandResolver.getSingleValue(args[0], ec.getRowIndex(), ec.getColumnIndex()) ;
			final ValueEval mu = OperandResolver.getSingleValue(args[1], ec.getRowIndex(), ec.getColumnIndex()) ;
			final ValueEval sigma = OperandResolver.getSingleValue(args[2], ec.getRowIndex(), ec.getColumnIndex()) ;			
			return new NumberEval(this.compute(OperandResolver.coerceValueToDouble(p), OperandResolver.coerceValueToDouble(mu), OperandResolver.coerceValueToDouble(sigma)));
		} catch (EvaluationException e) {
			return e.getErrorEval();
		}			 		
	}
 
	/**
	 * Original C++ implementation found at http://www.wilmott.com/messageview.cfm?catid=10&threadid=38771
	 * C# implementation found at http://weblogs.asp.net/esanchez/archive/2010/07/29/a-quick-and-dirty-implementation-of-excel-norminv-function-in-c.aspx
	 *    Compute the quantile function for the normal distribution.
	 *
	 *    For small to moderate probabilities, algorithm referenced
	 *    below is used to obtain an initial approximation which is
	 *    polished with a final Newton step.
	 *    For very large arguments, an algorithm of Wichura is used.
	 *
	 *  REFERENCE
	 *
	 *    Beasley, J. D. and S. G. Springer (1977).
	 *    Algorithm AS 111: The percentage points of the normal distribution,
	 *    Applied Statistics, 26, 118-121.
	 *
	 *     Wichura, M.J. (1988).
	 *     Algorithm AS 241: The Percentage Points of the Normal Distribution.
	 *     Applied Statistics, 37, 477-484.
	 * @param p
	 * @param mu
	 * @param sigma
	 * @return
	 */
	public double compute(double p, double mu, double sigma) {
		if(p < 0 || p > 1) 
			throw new RuntimeException("The probality p must be bigger than 0 and smaller than 1");		
		if(sigma < 0)
			throw new RuntimeException("The standard deviation sigma must be positive");
		if(p == 0)
			return Double.NEGATIVE_INFINITY;		
		if(p == 1)
			return Double.POSITIVE_INFINITY;		
		if(sigma == 0)
			return mu;		
		double  q, r, val;
 
		q = p - 0.5;
 
		/* 0.075 <= p <= 0.925 */
		if(Math.abs(q) <= .425) {
			r = .180625 - q * q;
			val =
		         q * (((((((r * 2509.0809287301226727 +
		                    33430.575583588128105) * r + 67265.770927008700853) * r +
		                  45921.953931549871457) * r + 13731.693765509461125) * r +
		                1971.5909503065514427) * r + 133.14166789178437745) * r +
		              3.387132872796366608)
		         / (((((((r * 5226.495278852854561 +
		                  28729.085735721942674) * r + 39307.89580009271061) * r +
		                21213.794301586595867) * r + 5394.1960214247511077) * r +
		              687.1870074920579083) * r + 42.313330701600911252) * r + 1);
		}
		/* closer than 0.075 from {0,1} boundary */
		else {
		     /* r = min(p, 1-p) < 0.075 */
			 if (q > 0) {
				 r = 1 - p;
			 } else {
				 r = p;
			 }
 
			 r = Math.sqrt(-Math.log(r));
			 /* r = sqrt(-log(r))  < ==>  min(p, 1-p) = exp( - r^2 ) */
 
			 if (r < = 5) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
				 r += -1.6;
				 val = (((((((r * 7.7454501427834140764e-4 +
		                     .0227238449892691845833) * r + .24178072517745061177) *
		                   r + 1.27045825245236838258) * r +
		                  3.64784832476320460504) * r + 5.7694972214606914055) *
		                r + 4.6303378461565452959) * r +
		               1.42343711074968357734)
		              / (((((((r *
		                       1.05075007164441684324e-9 + 5.475938084995344946e-4) *
		                      r + .0151986665636164571966) * r +
		                     .14810397642748007459) * r + .68976733498510000455) *
		                   r + 1.6763848301838038494) * r +
		                  2.05319162663775882187) * r + 1);
			 } else { /* very close to  0 or 1 */
				 r += -5;
				 val = (((((((r * 2.01033439929228813265e-7 +
		                     2.71155556874348757815e-5) * r +
		                    .0012426609473880784386) * r + .026532189526576123093) *
		                  r + .29656057182850489123) * r +
		                 1.7848265399172913358) * r + 5.4637849111641143699) *
		               r + 6.6579046435011037772)
		              / (((((((r *
		                       2.04426310338993978564e-15 + 1.4215117583164458887e-7) *
		                      r + 1.8463183175100546818e-5) * r +
		                     7.868691311456132591e-4) * r + .0148753612908506148525)
		                   * r + .13692988092273580531) * r +
		                  .59983220655588793769) * r + 1);
		      }
 
		      if (q < 0.0) {
		          val = -val;
		      }
		  }
 
		  return mu + sigma * val;		
	}
}

import org.apache.poi.ss.formula.OperationEvaluationContext; import org.apache.poi.ss.formula.eval.EvaluationException; import org.apache.poi.ss.formula.eval.NumberEval; import org.apache.poi.ss.formula.eval.OperandResolver; import org.apache.poi.ss.formula.eval.ValueEval; import org.apache.poi.ss.formula.functions.FreeRefFunction; /** * Java Implementation of http://support.microsoft.com/kb/827358/en-us implementing * Apache POIs {@link FreeRefFunction} to be used as a user defined function. * * Mathematics found here: * https://gist.github.com/kmpm/1211922/ * https://gist.github.com/kmpm/1211922/raw/a11e0dfc9fab493bcdadc669f3213d11f1897ebf/norminv.js * * Register for a given workbook with: * <code> final UDFFinder udfs = new DefaultUDFFinder(new String[]{ "MS_NormInv" }, new FreeRefFunction[]{ new NormInv() }) ; workbook.addToolPack(new AggregatingUDFFinder(new UDFFinder[] {udfs})); </code> * @author michael.simons */ public class NormInv implements FreeRefFunction { public ValueEval evaluate(ValueEval[] args, OperationEvaluationContext ec) { try { final ValueEval p = OperandResolver.getSingleValue(args[0], ec.getRowIndex(), ec.getColumnIndex()) ; final ValueEval mu = OperandResolver.getSingleValue(args[1], ec.getRowIndex(), ec.getColumnIndex()) ; final ValueEval sigma = OperandResolver.getSingleValue(args[2], ec.getRowIndex(), ec.getColumnIndex()) ; return new NumberEval(this.compute(OperandResolver.coerceValueToDouble(p), OperandResolver.coerceValueToDouble(mu), OperandResolver.coerceValueToDouble(sigma))); } catch (EvaluationException e) { return e.getErrorEval(); } } /** * Original C++ implementation found at http://www.wilmott.com/messageview.cfm?catid=10&threadid=38771 * C# implementation found at http://weblogs.asp.net/esanchez/archive/2010/07/29/a-quick-and-dirty-implementation-of-excel-norminv-function-in-c.aspx * Compute the quantile function for the normal distribution. * * For small to moderate probabilities, algorithm referenced * below is used to obtain an initial approximation which is * polished with a final Newton step. * For very large arguments, an algorithm of Wichura is used. * * REFERENCE * * Beasley, J. D. and S. G. Springer (1977). * Algorithm AS 111: The percentage points of the normal distribution, * Applied Statistics, 26, 118-121. * * Wichura, M.J. (1988). * Algorithm AS 241: The Percentage Points of the Normal Distribution. * Applied Statistics, 37, 477-484. * @param p * @param mu * @param sigma * @return */ public double compute(double p, double mu, double sigma) { if(p < 0 || p > 1) throw new RuntimeException("The probality p must be bigger than 0 and smaller than 1"); if(sigma < 0) throw new RuntimeException("The standard deviation sigma must be positive"); if(p == 0) return Double.NEGATIVE_INFINITY; if(p == 1) return Double.POSITIVE_INFINITY; if(sigma == 0) return mu; double q, r, val; q = p - 0.5; /* 0.075 <= p <= 0.925 */ if(Math.abs(q) <= .425) { r = .180625 - q * q; val = q * (((((((r * 2509.0809287301226727 + 33430.575583588128105) * r + 67265.770927008700853) * r + 45921.953931549871457) * r + 13731.693765509461125) * r + 1971.5909503065514427) * r + 133.14166789178437745) * r + 3.387132872796366608) / (((((((r * 5226.495278852854561 + 28729.085735721942674) * r + 39307.89580009271061) * r + 21213.794301586595867) * r + 5394.1960214247511077) * r + 687.1870074920579083) * r + 42.313330701600911252) * r + 1); } /* closer than 0.075 from {0,1} boundary */ else { /* r = min(p, 1-p) < 0.075 */ if (q > 0) { r = 1 - p; } else { r = p; } r = Math.sqrt(-Math.log(r)); /* r = sqrt(-log(r)) < ==> min(p, 1-p) = exp( - r^2 ) */ if (r < = 5) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */ r += -1.6; val = (((((((r * 7.7454501427834140764e-4 + .0227238449892691845833) * r + .24178072517745061177) * r + 1.27045825245236838258) * r + 3.64784832476320460504) * r + 5.7694972214606914055) * r + 4.6303378461565452959) * r + 1.42343711074968357734) / (((((((r * 1.05075007164441684324e-9 + 5.475938084995344946e-4) * r + .0151986665636164571966) * r + .14810397642748007459) * r + .68976733498510000455) * r + 1.6763848301838038494) * r + 2.05319162663775882187) * r + 1); } else { /* very close to 0 or 1 */ r += -5; val = (((((((r * 2.01033439929228813265e-7 + 2.71155556874348757815e-5) * r + .0012426609473880784386) * r + .026532189526576123093) * r + .29656057182850489123) * r + 1.7848265399172913358) * r + 5.4637849111641143699) * r + 6.6579046435011037772) / (((((((r * 2.04426310338993978564e-15 + 1.4215117583164458887e-7) * r + 1.8463183175100546818e-5) * r + 7.868691311456132591e-4) * r + .0148753612908506148525) * r + .13692988092273580531) * r + .59983220655588793769) * r + 1); } if (q < 0.0) { val = -val; } } return mu + sigma * val; } }

Even more bonus languages, here for Oracle PL/SQL:

CREATE OR REPLACE FUNCTION f_norminv(
	p     IN NUMBER,
	mu    IN NUMBER, 
	sigma IN NUMBER
) RETURN NUMBER
IS
	q   NUMBER;
	r   NUMBER;
	val NUMBER;
BEGIN
	IF p < 0 OR p > 1 THEN
		raise_application_error(-20100, 'The probality p must be bigger than 0 and smaller than 1');
	END IF;
	IF sigma < 0 THEN
		raise_application_error(-20100, 'The standard deviation sigma must be positive');
	END IF;
	IF p = 0 THEN 
		RETURN to_binary_double('-INF');	
	END IF;
	IF p = 1 THEN
		RETURN to_binary_double('INF');
	END IF;
	IF sigma = 0 THEN
		RETURN mu;
	END IF;
 
	q := p - 0.5;
 
	IF(ABS(q) <= .425) THEN
			r := .180625 - q * q;
			val :=
		         q * (((((((r * 2509.0809287301226727 +
		                    33430.575583588128105) * r + 67265.770927008700853) * r +
		                  45921.953931549871457) * r + 13731.693765509461125) * r +
		                1971.5909503065514427) * r + 133.14166789178437745) * r +
		              3.387132872796366608)
		         / (((((((r * 5226.495278852854561 +
		                  28729.085735721942674) * r + 39307.89580009271061) * r +
		                21213.794301586595867) * r + 5394.1960214247511077) * r +
		              687.1870074920579083) * r + 42.313330701600911252) * r + 1);
	ELSE
		/* r = min(p, 1-p) < 0.075 */
		IF q > 0 THEN
			r := 1 - p;
		ELSE
			r := p;
		END IF;
 
		r := SQRT(-LN(r));
		/* r = sqrt(-log(r))  <==>  min(p, 1-p) = exp( - r^2 ) */
 
		IF (r <= 5) THEN
			r := r - 1.6;
			val := (((((((r * 7.7454501427834140764e-4 +
		                 .0227238449892691845833) * r + .24178072517745061177) *
		               r + 1.27045825245236838258) * r +
		              3.64784832476320460504) * r + 5.7694972214606914055) *
		            r + 4.6303378461565452959) * r +
		           1.42343711074968357734)
		          / (((((((r *
		                   1.05075007164441684324e-9 + 5.475938084995344946e-4) *
		                  r + .0151986665636164571966) * r +
		                 .14810397642748007459) * r + .68976733498510000455) *
		               r + 1.6763848301838038494) * r +
		              2.05319162663775882187) * r + 1);
		ELSE /* very close to  0 or 1 */
			r := r - 5;
			val := (((((((r * 2.01033439929228813265e-7 +
		                 2.71155556874348757815e-5) * r +
		                .0012426609473880784386) * r + .026532189526576123093) *
		              r + .29656057182850489123) * r +
		             1.7848265399172913358) * r + 5.4637849111641143699) *
		           r + 6.6579046435011037772)
		          / (((((((r *
		                   2.04426310338993978564e-15 + 1.4215117583164458887e-7) *
		                  r + 1.8463183175100546818e-5) * r +
		                 7.868691311456132591e-4) * r + .0148753612908506148525)
		               * r + .13692988092273580531) * r +
		              .59983220655588793769) * r + 1);
		END IF;
 
		IF q < 0.0 THEN
			val := -val;
		END IF;
	END IF;
 
	RETURN mu + sigma * val;
END;
/
 
sho err;

create or replace function f_norminv( p IN NUMBER, mu IN NUMBER, sigma IN NUMBER ) RETURN NUMBER IS q NUMBER; r NUMBER; val NUMBER; BEGIN if p < 0 or p > 1 then raise_application_error(-20100, 'The probality p must be bigger than 0 and smaller than 1'); end if; if sigma < 0 then raise_application_error(-20100, 'The standard deviation sigma must be positive'); end if; if p = 0 then return to_binary_double('-INF'); end if; if p = 1 then return to_binary_double('INF'); end if; if sigma = 0 then return mu; end if; q := p - 0.5; if(abs(q) <= .425) then r := .180625 - q * q; val := q * (((((((r * 2509.0809287301226727 + 33430.575583588128105) * r + 67265.770927008700853) * r + 45921.953931549871457) * r + 13731.693765509461125) * r + 1971.5909503065514427) * r + 133.14166789178437745) * r + 3.387132872796366608) / (((((((r * 5226.495278852854561 + 28729.085735721942674) * r + 39307.89580009271061) * r + 21213.794301586595867) * r + 5394.1960214247511077) * r + 687.1870074920579083) * r + 42.313330701600911252) * r + 1); else /* r = min(p, 1-p) < 0.075 */ if q > 0 then r := 1 - p; else r := p; end if; r := sqrt(-ln(r)); /* r = sqrt(-log(r)) <==> min(p, 1-p) = exp( - r^2 ) */ if (r <= 5) then r := r - 1.6; val := (((((((r * 7.7454501427834140764e-4 + .0227238449892691845833) * r + .24178072517745061177) * r + 1.27045825245236838258) * r + 3.64784832476320460504) * r + 5.7694972214606914055) * r + 4.6303378461565452959) * r + 1.42343711074968357734) / (((((((r * 1.05075007164441684324e-9 + 5.475938084995344946e-4) * r + .0151986665636164571966) * r + .14810397642748007459) * r + .68976733498510000455) * r + 1.6763848301838038494) * r + 2.05319162663775882187) * r + 1); else /* very close to 0 or 1 */ r := r - 5; val := (((((((r * 2.01033439929228813265e-7 + 2.71155556874348757815e-5) * r + .0012426609473880784386) * r + .026532189526576123093) * r + .29656057182850489123) * r + 1.7848265399172913358) * r + 5.4637849111641143699) * r + 6.6579046435011037772) / (((((((r * 2.04426310338993978564e-15 + 1.4215117583164458887e-7) * r + 1.8463183175100546818e-5) * r + 7.868691311456132591e-4) * r + .0148753612908506148525) * r + .13692988092273580531) * r + .59983220655588793769) * r + 1); end if; if q < 0.0 then val := -val; end if; end if; return mu + sigma * val; END; / sho err;

7 comments

Post a comment | Trackback URI | RSS 2.0 feed for those comments

Yevgen wrote:

There was a mistake during translation to PL/SQL
After line “r := p;” should be “END IF;”
And last “END IF;” in function should be removed.

Posted on December 9, 2013 at 7:49 PM | Permalink
Michael wrote:

Hi Yevgen,
thank you very much, you’re absolutely right.

Post is updated.

Posted on December 9, 2013 at 9:42 PM | Permalink
Sid wrote:

Michael,
Your post is very useful.Thanks for providing PLSQL as well. do you have Java/PLSQL implementation of NORMSINV.?

Thanks in advance.!

Posted on September 28, 2016 at 3:02 PM | Permalink
Michael wrote:

Hey Sid,
yeah, it’s right there in the post… 🙂

Extract the #compute Method from the Java Class above into a a class that doesn’t depend on HSSF, something like
public class NormInv { public static double compute(double p, double mu, double sigma) { return 23.42; } }
public class NormInv { public static double compute(double p, double mu, double sigma) { return 23.42; } }
and follow the guide here https://docs.oracle.com/cd/B19306_01/java.102/b14187/chfive.htm#BABGEJDI

Posted on September 28, 2016 at 3:07 PM | Permalink
Sid wrote:

Thanks for Reply…i will work around that and come back in case need your help.!!

Posted on September 29, 2016 at 8:33 AM | Permalink
chinmaya wrote:

Thanks Michael . Your post is very very useful.

Posted on June 5, 2017 at 12:16 PM | Permalink
ömer Doğan wrote:

why i can’t implement android studio java

Posted on February 14, 2022 at 1:52 PM | Permalink

info.michael-simons.eu

Java implementation of Excels statistical functions NORMINV

7 comments

Post a Comment

Meta

Related posts

More from my site