Lazarus

Programming => Packages and Libraries => FPSpreadsheet => Topic started by: howardpc on May 10, 2021, 08:57:43 pm

Title: NORMSINV Excel function emulation
Post by: howardpc 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 (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.
Title: Re: NORMSINV Excel function emulation
Post by: Nitorami 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
Title: Re: NORMSINV Excel function emulation
Post by: wp 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.

Title: Re: NORMSINV Excel function emulation
Post by: howardpc 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.
Title: Re: NORMSINV Excel function emulation
Post by: wp 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