Normal Distributions in Delphi - Part II

This issue we continue our look at developing Statistical Routines to use in your Delphi applications, and in particular conclude looking at using Normal distributions.  These will be designed to use open array parameters where possible so that you can use them for standard arrays or with the dynamic arrays of Delphi 4 and later. While our freeware ESBMaths unit contains many stats routines (and more being added), if you want a good commercial stats library for Delphi, I would recommend either our ESBPCS or Turbo Power's SysTools 3.

In a latter issue, we will look at Math and Stats resources available on the Internet.

Reference

Formulae and other information used in this article rely on Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables by Milton Abramowitz and Irene A. Stegun, Dover. ISBN 0-486-61272-4. While not for the faint hearted this is a worthwhile book for Programmers interested in Mathematics.

Another Problem

After much research, it is found that the average life of the GloBrite Light Bulb under continuous use is 75 days, with a standard deviation of 15 days. GloBrite offers a free replacement if a Light Bulb expires within 50 days of purchase. What percentage of Light Bulbs would expire resulting in a free replacement?

Given the routines we encountered in Part 1 we could easily write code to handle the above:

    Value := NormalDistP (75, 15, 50);

Which would tell us that about 4.78% of the light bulbs sold would get a free replacement.

However, what if decreased production costs now mean that GloBrite is prepared to replace 7.5% of all bulbs sold under their replacement policy - to what could they extend their period for free replacement?

In other words, we want to be able to do the reverse of the above.

Inverse Normal Distribution

When it comes to probability distributions we normally need routines to calculate the Probabilities as we did last issue and we also need routines to calculate the Inverses. Most Approximations for Inverses of Distributions are only accurate to a few decimal places, though this is often more than sufficient.

As mentioned last issue, Abramowitz & Stegun provides many different Approximations that can be used and in particular the following is one for the Inverse Normal Distribution adapted into Delphi:

function InvNormalQ (const P: Extended): Extended;
const
  C0: Extended = 2.515517;
  C1: Extended = 0.802853;
  C2: Extended = 0.010328;
  D1: Extended = 1.432788;
  D2: Extended = 0.189269;
  D3: Extended = 0.001308;
var
  T: Extended;
begin
  T := Sqrt (Ln (1.0 / Sqr (P)));
  Result := T - (C0 + C1 * T + C2 * Sqr (T)) /
    (1.0 + D1 * T + D2 * Sqr (T) + D3 * Sqr (T) * T);
end;

The above works for probabilities that are greater than 0 and less than or equal to 0.5. The resultant error is less than 0.00045 - thus giving us 3 decimal place accuracy.

We now add a routine with a more useful interface and Exception Handling: 

{ Returns the Original AVal for Pr (X < AVal) = PVal if Less is
    True, or Pr (X > AVal) = PVal if Less is False,
    for a Normal Distribution with given Mean and Standard Deviation.

    Standard Deviation must be > 0 or function will return an Exception.
    PVal (Probability) must be 0 < PVal < 1 or function will return an Exception.

    Approximation has an absolute error < 4.5 x 10^4
    Thus 3 figure accuracy.

}
function InvNormalDist (const Mean, StdDev, PVal: Extended; const Less: Boolean): Extended;
var
  P: Extended;
  Lower: Boolean;
begin
  if not FloatIsPositive (StdDev) then
    raise EMathError.Create ('Standard Deviation must be Positive');

  if (PVal <= 0) or (PVal >= 1) then
    raise EMathError.Create ('Probability must be strictly between 0 and 1');

  if Less then
  begin
    P := 1 - PVal;
    Lower := P > 0.5;
    if Lower then
      P := PVal;
  end
  else
  begin
    P := PVal;
    Lower := P > 0.5;
    if Lower then
      P := 1 - PVal;
  end;
  Result := InvNormalQ (P);

  if Lower then
    Result := Mean - (Result * StdDev)
  else
    Result := Mean + (Result * StdDev);
end;

Thus we can now answer the problem that was proposed:

    Value := InvNormalDist (75, 15, 0.75, True);

Which results in finding out that the replacement period could be extended to 53 days.

MathRtns and Demo

From this Issue we now have a Zip file called MathsCorner.Zip - that will be updated each issue. It contains a Unit called MathRtns.Pas which contains all the routines so far covered in the Maths Corner as well as some extras.

You will also find NormDist2 Project which demonstrates the above problem.

Conclusion

Next Issue we will look at Angles.