Recent

Author Topic: NORMSINV Excel function emulation  (Read 1418 times)

howardpc

  • Hero Member
  • *****
  • Posts: 3753
NORMSINV Excel function emulation
« on: May 10, 2021, 08:57:43 pm »
I made a stab at porting  some C code I found on the internet for a NORMSINV function (which fpspreadsheet currently lacks) to Pascal.
The C code was based on an algorithm by Peter J Acklam (when he was associated with Oslo University).
Unfortunately the code blows up at the Ln() call. The quoted referenceURL:         http://www.math.uio.no/~jacklam
no longer exists, so I have not been able to consult the original algorithm to see where the code implementation is wrong.
Either the C code is faulty, or my port of it is wrong.
Can someone who understands normal distribution statistical functions point out where the error(s) are in the following?
Code: Pascal  [Select][+][-]
  1. uses  Math;
  2. function NormSInv(aValue, aMu, aSigma: Double): Double; // aMu = mean, aSigma = StdDev
  3. const
  4.   a: array[1..6] of Double = (-39.6968302866538,    220.946098424521,   -275.928510446969,
  5.                                138.357751867269,    -30.6647980661472,  2.50662827745924);
  6.   b: array[1..5] of Double = (-54.4760987982241,    161.585836858041,   -155.698979859887,
  7.                                66.8013118877197,    -13.2806815528857);
  8.   c: array[1..6] of Double = (-0.00778489400243029, -0.322396458041136, -2.40075827716184,
  9.                               -2.54973253934373,    4.37466414146497,   2.93816398269878);
  10.   d: array[1..4] of Double = (0.007784695709041462, 0.3224671290700398,
  11.                               2.445134137142996,    3.754408661907416);
  12.   pLow = 0.02425;
  13. var
  14.   q, r: Double;
  15. begin
  16.   p := (aValue - aMu)/aSigma;
  17.   Assert((p > 0.0) and (p < 1.0), 'NormDist: invalid first parameter');
  18.   q := Min(p, 1-p);
  19.   if q > pLow then
  20.     begin
  21.       q := q - 0.5;
  22.       r := q * q;
  23.       Result := q*(((((a[1] * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * r + a[6])/
  24.                   (((((b[1] * r + b[2]) * r + b[3]) * r + b[4]) * r + b[5]) * r + 1);
  25.     end
  26.   else
  27.     begin
  28.       q := Sqrt(-2 * Ln(q));
  29.       Result := (((((c[1] * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) * q + c[6]) /
  30.                  ((((d[1] * q + d[2]) * q + d[3]) * q + d[4]) * q + 1);
  31.     end
  32. end;                              
  33.  

Nitorami

  • Sr. Member
  • ****
  • Posts: 410
Re: NORMSINV Excel function emulation
« Reply #1 on: May 10, 2021, 09:08:45 pm »
No, but I have ported a version a while ago which works. It is part of a Class but should work if you simply remove the class.

Code: Pascal  [Select][+][-]
  1. class function StdNormal.iCDF (x: double): double;
  2. //Sehr exakte Version von Peter Acklam
  3. const
  4.   a1 = -3.969683028665376e+01;
  5.   a2 =  2.209460984245205e+02;
  6.   a3 = -2.759285104469687e+02;
  7.   a4 =  1.383577518672690e+02;
  8.   a5 = -3.066479806614716e+01;
  9.   a6 =  2.506628277459239e+00;
  10.  
  11.   b1 = -5.447609879822406e+01;
  12.   b2 =  1.615858368580409e+02;
  13.   b3 = -1.556989798598866e+02;
  14.   b4 =  6.680131188771972e+01;
  15.   b5 = -1.328068155288572e+01;
  16.  
  17.   c1 = -7.784894002430293e-03;
  18.   c2 = -3.223964580411365e-01;
  19.   c3 = -2.400758277161838e+00;
  20.   c4 = -2.549732539343734e+00;
  21.   c5 =  4.374664141464968e+00;
  22.   c6 =  2.938163982698783e+00;
  23.  
  24.   d1 =  7.784695709041462e-03;
  25.   d2 =  3.224671290700398e-01;
  26.   d3 =  2.445134137142996e+00;
  27.   d4 =  3.754408661907416e+00;
  28.  
  29. //Define break-points.
  30.   p_low  = 0.02425;
  31.   p_high = 1 - p_low;
  32.   eps = 1e-13;
  33. var
  34.   y,q,r: double;
  35.  
  36. //Rational approximation for lower region.
  37. //Selbe Formel wie upper region nur x transformiert und y negiert.
  38. begin
  39.   assert ((x>0) and (x<1),'stdnormal.icdf: invalid argument outside ]0..1[');
  40.   if x < p_low then begin
  41.      q := sqrt(-2*ln(x));
  42.      y := (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
  43.   end
  44.   //Rational approximation for central region.
  45.   else if x <= p_high then begin
  46.      q := x - 0.5;
  47.      r := q*q;
  48.      y := (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1);
  49.   end
  50.   //Rational approximation for upper region.
  51.   else begin
  52.      q := sqrt(-2*ln(1-x));
  53.      y := -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
  54.   end;
  55.   result := y;
  56. end;
  57.  

edit: icdf (x,mu,sigma) = icdf(x) *sigma + mu
« Last Edit: May 10, 2021, 10:02:08 pm by Nitorami »

wp

  • Hero Member
  • *****
  • Posts: 8727
Re: NORMSINV Excel function emulation
« Reply #2 on: May 10, 2021, 10:30:45 pm »
Howard, I cannot help you, I cannot find that reference even in the Wayback Machine (*). But if you need the function you could also use the invnormaldist() in the "spe" unit of fpc's NumLib.

I wanted to keep NumLib out of fpspreadsheet, and therefore, I do not plan to add special funtions. But the fpspreadsheet formula system is open and can be extended by any user formula - see the example in folder examples/other/user_defined_formulas.

----------------
(*) Found this: https://stackedboxes.org/2017/05/01/acklams-normal-quantile-function/, which seems to be Nitorami's code.

« Last Edit: May 10, 2021, 10:35:50 pm by wp »
Mainly Lazarus trunk / fpc 3.2.0 / all 32-bit on Win-10, but many more...

howardpc

  • Hero Member
  • *****
  • Posts: 3753
Re: NORMSINV Excel function emulation
« Reply #3 on: May 10, 2021, 10:42:45 pm »
Thanks to both of you. Nitorami's code works well.
I did not know about the spe unit, so I'm glad you pointed it out.

wp

  • Hero Member
  • *****
  • Posts: 8727
Re: NORMSINV Excel function emulation
« Reply #4 on: May 10, 2021, 11:04:15 pm »
I did not know about the spe unit, so I'm glad you pointed it out.
Then maybe I should also point you to the wiki help of NumLib: https://wiki.lazarus.freepascal.org/NumLib_Documentation
Mainly Lazarus trunk / fpc 3.2.0 / all 32-bit on Win-10, but many more...

 

TinyPortal © 2005-2018