TAChart: Add TCumulNormDistrAxisTransform based on code by Werner Pamler. Add test.

git-svn-id: trunk@32691 -

components/tachart/tamath.pas

@@ -11,7 +11,7 @@
  *                                                                           *
  *****************************************************************************
 
- Authors: Alexander Klenin
+ Authors: Alexander Klenin, Werner Pamler
 
 }
 unit TAMath;
@@ -23,6 +23,9 @@ interface
 uses
   Classes, SysUtils; 
 
+function CumulNormDistr(AX: Double): Double;
+function InvCumulNormDistr(AX: Double): Double;
+
 procedure EnsureOrder(var A, B: Integer); overload; inline;
 procedure EnsureOrder(var A, B: Double); overload; inline;
 
@@ -34,7 +37,87 @@ function SafeNan: Double; inline;
 implementation
 
 uses
-  Math, TAChartUtils;
+  Math, spe, TAChartUtils;
+
+// Cumulative normal distribution
+// x = -INF ... INF --> Result = 0 ... 1
+function CumulNormDistr(AX: Double): Double;
+begin
+  if AX > 0 then
+    Result := (speerf(AX / Sqrt(2)) + 1) * 0.5
+  else if AX < 0 then
+    Result := (1 - speerf(-AX / Sqrt(2))) * 0.5
+  else
+    Result := 0;
+end;
+
+// Inverse cumulative normal distribution.
+// x = 0 ... 1 --> Result = -INF ... +INF
+// Algorithm by Peter John Acklam.
+// http://home.online.no/~pjacklam/notes/invnorm/index.html
+function InvCumulNormDistr(AX: Double): Double;
+const
+  A: array[1..6] of Double = (
+    -3.969683028665376e+01,
+    +2.209460984245205e+02,
+    -2.759285104469687e+02,
+    +1.383577518672690e+02,
+    -3.066479806614716e+01,
+    +2.506628277459239e+00
+  );
+  B: array[1..5] of Double = (
+    -5.447609879822406e+01,
+    +1.615858368580409e+02,
+    -1.556989798598866e+02,
+    +6.680131188771972e+01,
+    -1.328068155288572e+01
+  );
+  C: array[1..6] of Double = (
+    -7.784894002430293e-03,
+    -3.223964580411365e-01,
+    -2.400758277161838e+00,
+    -2.549732539343734e+00,
+    +4.374664141464968e+00,
+    +2.938163982698783e+00
+  );
+  D: array[1..4] of Double = (
+    +7.784695709041462e-03,
+    +3.224671290700398e-01,
+    +2.445134137142996e+00,
+    +3.754408661907416e+00
+  );
+  // Switching points between regions.
+  P_LOW = 0.02425;
+  P_HIGH = 1 - P_LOW;
+var
+  q, r: Extended;
+begin
+  if AX <= 0 then
+    Result := NegInfinity
+  else if AX < P_LOW then begin
+    // Rational approximation for lower region.
+    q := Sqrt(-2 * Ln(AX));
+    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
+  else if AX <= P_HIGH then begin
+    // Rational approximation for central region.
+    q := AX - 0.5 ;
+    r := q * q ;
+    Result :=
+      (((((A[1] * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * r + A[6]) * q /
+      (((((B[1] * r + B[2]) * r + B[3]) * r + B[4]) * r + B[5]) * r + 1);
+  end
+  else if AX < 1 then begin
+    // Rational approximation for upper region.
+    q := Sqrt(-2 * Ln(1 - AX));
+    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 else
+    Result := SafeInfinity;
+end;
 
 procedure EnsureOrder(var A, B: Integer); overload; inline;
 begin

components/tachart/tatransformations.pas

@@ -33,7 +33,7 @@ type
   { TAxisTransform }
 
   TAxisTransform = class(TIndexedComponent)
-  private
+  strict private
     FEnabled: Boolean;
     FTransformations: TChartAxisTransformations;
     procedure SetEnabled(AValue: Boolean);
@@ -82,7 +82,7 @@ type
   { TChartAxisTransformations }
 
   TChartAxisTransformations = class(TComponent)
-  private
+  strict private
     FBroadcaster: TBroadcaster;
     FList: TAxisTransformList;
   protected
@@ -108,7 +108,7 @@ type
   { TLinearAxisTransform }
 
   TLinearAxisTransform = class(TAxisTransform)
-  private
+  strict private
     FOffset: Double;
     FScale: Double;
     function OffsetIsStored: Boolean;
@@ -130,7 +130,7 @@ type
   { TAutoScaleAxisTransform }
 
   TAutoScaleAxisTransform = class(TAxisTransform)
-  private
+  strict private
     FMaxValue: Double;
     FMinValue: Double;
     function MaxValueIsStored: Boolean;
@@ -158,7 +158,7 @@ type
   { TLogarithmAxisTransform }
 
   TLogarithmAxisTransform = class(TAxisTransform)
-  private
+  strict private
     FBase: Double;
     procedure SetBase(AValue: Double);
   public
@@ -172,6 +172,12 @@ type
     property Base: Double read FBase write SetBase;
   end;
 
+  TCumulNormDistrAxisTransform = class(TAxisTransform)
+  public
+    function AxisToGraph(AX: Double): Double; override;
+    function GraphToAxis(AX: Double): Double; override;
+  end;
+
   TTransformEvent = procedure (AX: Double; out AT: Double) of object;
 
   { TUserDefinedAxisTransform }
@@ -732,6 +738,18 @@ begin
     AMax := MaxValue;
 end;
 
+{ TCumulNormDistrAxisTransform }
+
+function TCumulNormDistrAxisTransform.AxisToGraph(AX: Double): Double;
+begin
+  Result := InvCumulNormDistr(AX);
+end;
+
+function TCumulNormDistrAxisTransform.GraphToAxis(AX: Double): Double;
+begin
+  Result := CumulNormDistr(AX);
+end;
+
 { TUserDefinedAxisTransform }
 
 procedure TUserDefinedAxisTransform.Assign(ASource: TPersistent);
@@ -778,6 +796,8 @@ initialization
 
   AxisTransformsClassRegistry := TStringList.Create;
   RegisterAxisTransformClass(TAutoScaleAxisTransform, 'Auto scale');
+  RegisterAxisTransformClass(
+    TCumulNormDistrAxisTransform, 'Cumulative normal distribution');
   RegisterAxisTransformClass(TLinearAxisTransform, 'Linear');
   RegisterAxisTransformClass(TLogarithmAxisTransform, 'Logarithmic');
   RegisterAxisTransformClass(TUserDefinedAxisTransform, 'User defined');

components/tachart/test/UtilsTest.pas

@@ -37,6 +37,10 @@ type
     procedure Merge;
   end;
 
+  TMathTest = class(TTestCase)
+  published
+    procedure CumulNurmDistrTest;
+  end;
 
   TGeometryTest = class(TTestCase)
   strict private
@@ -78,7 +82,7 @@ type
 implementation
 
 uses
-  TAGeometry;
+  TAGeometry, TAMath;
 
 { TIntervalListTest }
 
@@ -147,6 +151,20 @@ begin
   FreeAndNil(FIList);
 end;
 
+{ TMathTest }
+
+procedure TMathTest.CumulNurmDistrTest;
+const
+  INV_PTS: array [1..3] of Double = (-1.5, 0.33, 2.0);
+var
+  p: Double;
+begin
+  AssertEquals(0, CumulNormDistr(0));
+  AssertEquals(0.84134, CumulNormDistr(1.0));
+  for p in INV_PTS do
+    AssertEquals(p, InvCumulNormDistr(CumulNormDistr(p)));
+end;
+
 { TGeometryTest }
 
 procedure TGeometryTest.AssertEquals(const Expected, Actual: TDoublePoint);
@@ -414,7 +432,7 @@ end;
 initialization
 
   RegisterTests([
-    TIntervalListTest, TGeometryTest, TColorTest, TRTTITest,
+    TIntervalListTest, TMathTest, TGeometryTest, TColorTest, TRTTITest,
     TPublishedIntegerSetTest]);
 
 end.