NORMSINV Excel function emulation

(1/1)

howardpc:
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  [+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i < x.length; i++) { x[i].style.maxHeight='none'; x[i].style.height = Math.min(x[i].clientHeight+15,306)+'px'; x[i].style.resize = "vertical";}};} ---uses  Math;function NormSInv(aValue, aMu, aSigma: Double): Double; // aMu = mean, aSigma = StdDevconst  a: array[1..6] of Double = (-39.6968302866538,    220.946098424521,   -275.928510446969,                               138.357751867269,    -30.6647980661472,  2.50662827745924);  b: array[1..5] of Double = (-54.4760987982241,    161.585836858041,   -155.698979859887,                               66.8013118877197,    -13.2806815528857);  c: array[1..6] of Double = (-0.00778489400243029, -0.322396458041136, -2.40075827716184,                              -2.54973253934373,    4.37466414146497,   2.93816398269878);  d: array[1..4] of Double = (0.007784695709041462, 0.3224671290700398,                              2.445134137142996,    3.754408661907416);  pLow = 0.02425;var  q, r: Double;begin  p := (aValue - aMu)/aSigma;  Assert((p > 0.0) and (p < 1.0), 'NormDist: invalid first parameter');  q := Min(p, 1-p);  if q > pLow then    begin      q := q - 0.5;      r := q * q;      Result := q*(((((a[1] * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * r + a[6])/                  (((((b[1] * r + b[2]) * r + b[3]) * r + b[4]) * r + b[5]) * r + 1);    end  else    begin      q := Sqrt(-2 * Ln(q));      Result := (((((c[1] * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) * q + c[6]) /                 ((((d[1] * q + d[2]) * q + d[3]) * q + d[4]) * q + 1);    endend;

Nitorami:
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  [+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i < x.length; i++) { x[i].style.maxHeight='none'; x[i].style.height = Math.min(x[i].clientHeight+15,306)+'px'; x[i].style.resize = "vertical";}};} ---class function StdNormal.iCDF (x: double): double;//Sehr exakte Version von Peter Acklamconst  a1 = -3.969683028665376e+01;  a2 =  2.209460984245205e+02;  a3 = -2.759285104469687e+02;  a4 =  1.383577518672690e+02;  a5 = -3.066479806614716e+01;  a6 =  2.506628277459239e+00;   b1 = -5.447609879822406e+01;  b2 =  1.615858368580409e+02;  b3 = -1.556989798598866e+02;  b4 =  6.680131188771972e+01;  b5 = -1.328068155288572e+01;   c1 = -7.784894002430293e-03;  c2 = -3.223964580411365e-01;  c3 = -2.400758277161838e+00;  c4 = -2.549732539343734e+00;  c5 =  4.374664141464968e+00;  c6 =  2.938163982698783e+00;   d1 =  7.784695709041462e-03;  d2 =  3.224671290700398e-01;  d3 =  2.445134137142996e+00;  d4 =  3.754408661907416e+00; //Define break-points.  p_low  = 0.02425;  p_high = 1 - p_low;  eps = 1e-13;var  y,q,r: double; //Rational approximation for lower region.//Selbe Formel wie upper region nur x transformiert und y negiert.begin  assert ((x>0) and (x<1),'stdnormal.icdf: invalid argument outside ]0..1[');  if x < p_low then begin     q := sqrt(-2*ln(x));     y := (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);  end  //Rational approximation for central region.  else if x <= p_high then begin     q := x - 0.5;     r := q*q;     y := (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1);  end  //Rational approximation for upper region.  else begin     q := sqrt(-2*ln(1-x));     y := -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);  end;  result := y;end;
edit: icdf (x,mu,sigma) = icdf(x) *sigma + mu

wp:
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.

howardpc:
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:

--- Quote from: howardpc on May 10, 2021, 10:42:45 pm ---I did not know about the spe unit, so I'm glad you pointed it out.

--- End quote ---
Then maybe I should also point you to the wiki help of NumLib: https://wiki.lazarus.freepascal.org/NumLib_Documentation