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?
uses Math;
function NormSInv(aValue, aMu, aSigma: Double): Double; // aMu = mean, aSigma = StdDev
const
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);
end
end;