LazStats: Further refactoring of ResistanceLineUnit. Some rearrangement of units to make MathUnit only dependent on Globals (for easier testing).

git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7756 8e941d3f-bd1b-0410-a28a-d453659cc2b4

applications/lazstats/source/forms/analysis/descriptive/multxvsyunit.pas

@@ -67,7 +67,7 @@ implementation
 
 uses
   TATypes,
-  Math, Utils;
+  Math, Utils, MathUnit;
 
 
 { TMultXvsYFrm }

applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas

@@ -1,7 +1,7 @@
 // Use file "Sickness.laz" for testing
 
 // to do Not sure about the "red" slope line...
-// to do: There are also two positions where a negative intercept is made positive. Why? Correct?
+// to do: There are also two code positions where a negative intercept is made positive. Why? Correct?
 
 unit ResistanceLineUnit;
 
@@ -12,9 +12,8 @@ interface
 uses
   Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
   StdCtrls, ExtCtrls, Buttons, Printers, ComCtrls,
-  MainUnit, Globals, FunctionsLib, OutputUnit, DataProcs, DictionaryUnit,
+  MainUnit, Globals, DataProcs, DictionaryUnit,
-  //BlankFrmUnit,
+  MathUnit, BasicStatsReportAndChartFormUnit;
-  BasicStatsReportAndChartFormUnit;
 
 type
 
@@ -55,8 +54,8 @@ type
   private
     procedure PlotMedians(const XMedians, YMedians: DblDyneVec;
       ASlope, AIntercept: Double);
-    procedure PlotXY(const XPoints, YPoints, UpConf, LowConf: DblDyneVec;
+    procedure PlotXY(const XPoints, YPoints: DblDyneVec;
-      XMean, YMean, R, ASlope, AIntercept: Double);
+      ARegressionResults: TBivariateRegressionResults);
 
     procedure ResidualsToGrid(xCol, yCol: Integer;
       const AColNoSelected: IntDyneVec; Slope, Intercept: Double);
@@ -68,8 +67,7 @@ type
       const XMedians, YMedians: DblDyneVec; const GrpSize: IntDyneVec;
       Slope1, Slope2, Slope, Intercept: Double);
     procedure WriteXYReport(AReport: TStrings; AFileName: String;
-      xMean, yMean, xVariance, yVariance, xStdDev, yStdDev, R, Slope, Intercept, SEPred: double;
+      const ARegressionResults: TBivariateRegressionResults);
-      N: Integer);
 
   protected
     procedure AdjustConstraints; override;
@@ -91,7 +89,7 @@ implementation
 uses
   Math, Grids,
   TATypes, TAChartUtils, TAChartAxisUtils, TALegend, TASources, TACustomSeries,
-  ChartFrameUnit, Utils, MathUnit, GridProcs;
+  ChartFrameUnit, Utils, GridProcs;
 
 { TResistanceLineForm }
 
@@ -142,15 +140,13 @@ var
   xMedians: DblDyneVec = nil;
   yMedians: DblDyneVec = nil;
   colNoSelected : IntDyneVec = nil;
-  upConf: DblDyneVec = nil;
-  lowConf: DblDyneVec = nil;
   grpSize: IntDyneVec = nil;
   xCol, yCol: Integer;
-  i, N, DF: integer;
+  regressionRes: TBivariateRegressionResults;
-  xMean, yMean, xVariance, yVariance, xStdDev, yStdDev, SXX, SXY, SYY, t: Double;
+  N: Integer;
-  confBand, yPred: double;
+  confBand: Double;
-  R, SEPred, interceptFit, slopeFit, slope1, slope2, slopeRL, seData: double;
+  slope1, slope2, slopeRL: Double;
-  c, c1, c2, c3 : double; // constants obtained from control points
+  c1, c2, c3, c: Double;
   lReport: TStrings;
   dataGrid: TStringGrid;
 begin
@@ -178,9 +174,9 @@ begin
   xyPoints[0] := CollectValues(dataGrid, xCol, colNoSelected);
   xyPoints[1] := CollectValues(dataGrid, ycol, colNoSelected);
   N := Length(xyPoints[0]);
-  if N = 0 then
+  if N < 3 then
   begin
-    ErrorMsg('No data points.');
+    ErrorMsg('At least three data points required.');
     exit;
   end;
   if N <> Length(xyPoints[1]) then
@@ -192,36 +188,9 @@ begin
   // Sort on the x values
   SortOnX(xyPoints[0], xyPoints[1]);
 
-  // Calculate means, variances, stddevs
+  // Calculate bivariate regression
-  Calc_MeanVarStdDevSS(xyPoints[0], xMean, xVariance, xStdDev, SXX);
+  confBand := StrToFloat(ConfEdit.Text) / 100;
-  Calc_MeanVarStdDevSS(xyPoints[1], yMean, yVariance, yStdDev, SYY);
+  Calc_BivariateRegression(xyPoints[0], xyPoints[1], confBand, regressionRes);
-
-  // Calculate linear fit of x vs y
-  SXY := 0;
-  for i := 0 to N-1 do
-    SXY := SXY + xyPoints[0, i] * xyPoints[1, i];
-
-  R := (SXY - xMean * yMean * N) / ((N - 1) * xStdDev * yStdDev);
-  sePred := sqrt(1.0 - sqr(R)) * yStdDev * sqrt((N - 1) / (N - 2));
-  slopeFit := R * yStdDev / xStdDev;
-  interceptFit := yMean - slopeFit * xMean;
-
-  // Get upper and lower confidence limits for each x value
-  if ConfChk.Checked then
-  begin
-    SetLength(upConf, N);
-    SetLength(lowConf, N);
-    confBand := StrToFloat(ConfEdit.Text) / 100.0;
-    DF := N - 2;
-    t := InverseT(ConfBand, DF);
-    for i := 0 to N-1 do
-    begin
-      yPred := interceptFit + slopeFit * xyPoints[0, i];
-      seData := sePred * sqrt(1.0 + 1/N + sqr(xyPoints[0, i] - xMean)/SXX);
-      upConf[i] := yPred + t * seData;
-      lowConf[i] := yPred - t * seData;
-    end;
-  end;
 
   // Do the resistant line analysis
   ResistantLineAnalysis(xyPoints, xMedians, yMedians, grpSize);
@@ -240,6 +209,7 @@ begin
   if GridChk.Checked then
     ResidualsToGrid(xCol, yCol, colNoSelected, slopeRL, c);
 
+  // Write results to report in ReportFrame
   lReport := TStringList.Create;
   try
     WriteMedianReport(lReport,
@@ -250,9 +220,7 @@ begin
     if DescChk.Checked then
        WriteXYReport(lReport,
          OS3MainFrm.FileNameEdit.Text,
-         xMean, yMean, xVariance, yVariance, xStdDev, yStdDev,
+         regressionRes
-         R, slopeFit, interceptFit, sePred,
-         N
        );
 
     FReportFrame.DisplayReport(lReport);
@@ -260,8 +228,9 @@ begin
     lReport.Free;
   end;
 
+  // Plot results in ChartFrame
   if PlotPointsChk.Checked then
-    PlotXY(xyPoints[0], xyPoints[1], upConf, lowConf, xMean, yMean, R, slopeFit, interceptFit);
+    PlotXY(xyPoints[0], xyPoints[1], regressionRes);
 
   if PlotMediansChk.Checked then
     PlotMedians(XMedians, YMedians, slopeRL, c);
@@ -312,14 +281,16 @@ begin
 end;
 
 
-procedure TResistanceLineForm.PlotXY(const XPoints, YPoints, UpConf, LowConf: DblDyneVec;
+procedure TResistanceLineForm.PlotXY(const XPoints, YPoints: DblDyneVec;
-  XMean, YMean, R, ASlope, AIntercept: Double);
+  ARegressionResults: TBivariateRegressionResults);
 var
   ser: TChartSeries;
   xPts: DblDyneVec = nil;
   yPts: DblDyneVec = nil;
+  conf: DblDyneVec = nil;
   rightLabels: TListChartSource;
   topLabels: TListChartSource;
+  i: Integer;
 begin
   rightLabels := FChartFrame.Chart.AxisList[2].Marks.Source as TListChartSource;
   rightLabels.Clear;
@@ -331,19 +302,21 @@ begin
   FChartFrame.Chart.Margins.Right := 4;
 
   FChartFrame.SetTitle('X vs Y plot using file '+ OS3MainFrm.FileNameEdit.Text);
-  FChartFrame.SetFooter(Format('R(X,Y) = %.3f, Slope = %.2f, Intercept = %.2f', [R, ASlope, AIntercept]));
+  with ARegressionResults do
+    FChartFrame.SetFooter(Format('R(X,Y) = %.3f, Slope = %.2f, Intercept = %.2f', [R, Slope, Intercept]));
   FChartFrame.SetXTitle(XEdit.Text);
   FChartFrame.SetYTitle(YEdit.Text);
 
   // Draw means
   if MeansChk.Checked then
-  begin
+    with ARegressionResults do
-    FChartFrame.VertLine(XMean, clGreen, psDashDot, 'Mean ' + XEdit.Text);
+    begin
-    topLabels.Add(XMean, XMean, 'Mean ' + XEdit.Text);
+      FChartFrame.VertLine(XMean, clGreen, psDashDot, 'Mean ' + XEdit.Text);
+      topLabels.Add(XMean, XMean, 'Mean ' + XEdit.Text);
 
-    FChartFrame.HorLine(YMean, clGreen, psDash, 'Mean ' + YEdit.Text);
+      FChartFrame.HorLine(YMean, clGreen, psDash, 'Mean ' + YEdit.Text);
-    rightLabels.Add(YMean, YMean, 'Mean ' + YEdit.Text);
+      rightLabels.Add(YMean, YMean, 'Mean ' + YEdit.Text);
-  end;
+    end;
 
   // Draw slope line
   if LineChk.Checked then
@@ -351,24 +324,29 @@ begin
     SetLength(xPts, 2);
     SetLength(yPts, 2);
     xPts[0] := XPoints[0];
-    yPts[0] := AIntercept + XPoints[0] * ASlope;
     xPts[1] := XPoints[High(XPoints)];
-    yPts[1] := AIntercept + XPoints[High(XPoints)] * ASlope;
+    with ARegressionResults do
+    begin
+      yPts[0] := Intercept + XPoints[0] * Slope;
+      yPts[1] := Intercept + XPoints[High(XPoints)] * Slope;
+    end;
     ser := FChartFrame.PlotXY(ptLines, xpts, ypts, nil, nil, 'Regression', clBlack);
     rightLabels.Add(ser.XValue[ser.Count-1], ser.YValue[ser.Count-1], 'Regression');
   end;
 
-  // Draw upper confidence band
   if ConfChk.Checked then
   begin
-    ser := FChartFrame.PlotXY(ptLines, XPoints, UpConf, nil, nil, 'Upper confidence band', clRed);
+    // Draw upper confidence band
+    SetLength(conf, ARegressionResults.Count);
+    for i := 0 to ARegressionResults.Count-1 do
+      conf[i] := ARegressionResults.ConfidenceLimits(XPoints[i], true);
+    ser := FChartFrame.PlotXY(ptLines, XPoints, conf, nil, nil, 'Upper confidence band', clRed);
     rightLabels.Add(ser.yValue[ser.Count-1], ser.YValue[ser.Count-1], 'UCL');
-  end;
 
-  // Draw lower confidence band
+    // Draw lower confidence band
-  if ConfChk.Checked then
+    for i := 0 to ARegressionResults.Count-1 do
-  begin
+      conf[i] := ARegressionResults.ConfidenceLimits(XPoints[i], false);
-    ser := FChartFrame.PlotXY(ptLines, XPoints, LowConf, nil, nil, 'Lower confidence band', clRed);
+    ser := FChartFrame.PlotXY(ptLines, XPoints, conf, nil, nil, 'Lower confidence band', clRed);
     rightLabels.Add(ser.yValue[ser.Count-1], ser.YValue[ser.Count-1], 'LCL');
   end;
 
@@ -587,8 +565,11 @@ end;
 
 
 procedure TResistanceLineForm.WriteXYReport(AReport: TStrings; AFileName: String;
+  const ARegressionResults: TBivariateRegressionResults);
+{
   xMean, yMean, xVariance, yVariance, xStdDev, yStdDev, R, slope, intercept, SEPred: double;
   N: Integer);
+  }
 begin
   AReport.Add('');
   AReport.Add(DIVIDER);
@@ -599,10 +580,22 @@ begin
   AReport.Add('Variables:');
   AReport.Add('  X: ' + XEdit.Text);
   AReport.Add('  Y: ' + YEdit.Text);
-//  AReport.Add('X = %s, Y = %s from file: %s', [XEdit.Text, YEdit.Text, AFileName]);
   AReport.Add('');
   AReport.Add('Variable      Mean    Variance  Std.Dev.');
   AReport.Add('----------  --------  --------  --------');
+  with ARegressionResults do
+  begin
+    AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [XEdit.Text, XMean, XVariance, XStdDev]);
+    AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [YEdit.Text, YMean, YVariance, YStdDev]);
+    AReport.Add('');
+    AReport.Add('Correlation:                %8.4f', [R]);
+    AReport.Add('Slope:                      %8.3f', [Slope]);
+    AReport.Add('Intercept:                  %8.3f', [Intercept]);
+    AReport.Add('Standard Error of Estimate: %8.3f', [StdErrorPredicted]);
+    AReport.Add('Number of good cases:       %8d',   [Count]);
+  end;
+
+  {
   AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [XEdit.Text, xMean, xVariance, xStdDev]);
   AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [YEdit.Text, yMean, yVariance, yStdDev]);
   AReport.Add('');
@@ -611,6 +604,7 @@ begin
   AReport.Add('Intercept:                  %8.3f', [Intercept]);
   AReport.Add('Standard Error of Estimate: %8.3f', [SEPred]);
   AReport.Add('Number of good cases:       %8d', [N]);
+  }
 end;
 
 

applications/lazstats/source/forms/analysis/measurement_programs/gradebookunit.pas

@@ -105,7 +105,7 @@ implementation
 
 uses
   Math,
-  Utils;
+  Utils, MathUnit;
 
 { TGradebookFrm }
 

applications/lazstats/source/forms/analysis/measurement_programs/gradingunit.pas

@@ -72,7 +72,7 @@ implementation
 
 uses
   Math,
-  Utils, GradebookUnit;
+  Utils, MathUnit, GradebookUnit;
 
 { TGradingFrm }
 

applications/lazstats/source/forms/analysis/multivariate/kmeansunit.pas

@@ -85,7 +85,8 @@ var
 implementation
 
 uses
-  Math, Utils;
+  Math,
+  Utils, MathUnit;
 
 { TKMeansFrm }
 

applications/lazstats/source/forms/analysis/multivariate/medianpolishunit.pas

@@ -83,7 +83,8 @@ var
 implementation
 
 uses
-  Math, Utils;
+  Math,
+  Utils, MathUnit;
 
 { TMedianPolishForm }
 

applications/lazstats/source/forms/analysis/multivariate/pathunit.pas

@@ -84,7 +84,8 @@ var
 implementation
 
 uses
-  Math, Utils;
+  Math,
+  Utils, MathUnit;
 
 { TPathFrm }
 

applications/lazstats/source/forms/analysis/multivariate/singlelinkunit.pas

@@ -61,7 +61,8 @@ var
 implementation
 
 uses
-  Math, Utils;
+  Math,
+  Utils, MathUnit;
 
 { TSingleLinkFrm }
 

applications/lazstats/source/forms/analysis/nonparametric/sensunit.pas

@@ -65,7 +65,8 @@ var
 implementation
 
 uses
-  Math, Utils;
+  Math,
+  Utils, MathUnit;
 
 { TSensForm }
 

applications/lazstats/source/forms/misc/graphlib.pas

@@ -9,6 +9,15 @@ uses
   ExtCtrls, StdCtrls, Printers, Math,
   Globals;
 
+const
+  DATA_COLORS: array[0..10] of TColor = (
+    $A1A45D, $3153C4, $0996E7, $4AE8F6, $A7A2B1,
+    $84A7C9, $51798C, $87CDD8, $536508, $7BD8F3,
+    $846402
+  );
+  // NOTE: This is a duplication of the declaration in ChartUnitFrame and will
+  // become obsolete when the GraphLib unit will be finally replaced by TAChart.
+
 
 type
 

applications/lazstats/source/frames/chartframeunit.pas

@@ -10,6 +10,18 @@ uses
   TATools,
   Globals, MainDM;
 
+const
+  DATA_COLORS: array[0..10] of TColor = (
+    $A1A45D, $3153C4, $0996E7, $4AE8F6, $A7A2B1,
+    $84A7C9, $51798C, $87CDD8, $536508, $7BD8F3,
+    $846402
+  );
+  DATA_SYMBOLS: array[0..5] of TSeriesPointerStyle = (
+    psRectangle, psCircle, psDiamond,
+    psDownTriangle, psHexagon, psFullStar
+  );
+
+
 type
   TPlotType = (ptLines, ptSymbols, ptLinesAndSymbols, ptBars,
     ptArea);

applications/lazstats/source/units/dataprocs.pas

@@ -61,7 +61,7 @@ function StringsToInt(StrCol: integer; var NewCol: integer; Prompt: boolean): bo
 implementation
 
 uses
-  Utils, MainUnit;
+  Utils, MathUnit, MainUnit;
 
 
 // NOTE: Do not call GridProcs.GoodRecord here because this old function may

applications/lazstats/source/units/globals.pas

@@ -5,7 +5,7 @@ unit Globals;
 interface
 
 uses
-  Classes, SysUtils, Graphics, TATypes;
+  Classes, SysUtils;
 
 type
   IntDyneVec = array of integer;
@@ -75,27 +75,14 @@ const
   DEFAULT_ALPHA_LEVEL = 0.05;
   DEFAULT_BETA_LEVEL = 0.20;
 
-  {
-  DATA_COLORS: array[0..11] of TColor = (
-    clMaroon, clRed, clBlue, clGreen, clNavy, clTeal,
-    clAqua, clLime, clFuchsia, clGray, clSilver, clOlive
-  );
-  }
-  DATA_COLORS: array[0..10] of TColor = (
-    $A1A45D, $3153C4, $0996E7, $4AE8F6, $A7A2B1, $84A7C9, $51798C, $87CDD8, $536508, $7BD8F3, $846402
-  );
-  DATA_SYMBOLS: array[0..5] of TSeriesPointerStyle = (psRectangle, psCircle, psDiamond,
-     psDownTriangle, psHexagon, psFullStar);
-
-
   DIVIDER = '===========================================================================';
   DIVIDER_SMALL = '---------------------------------------------------------------------------';
   DIVIDER_AUTO = '@=';
   DIVIDER_SMALL_AUTO = '@-';
 
-  GRAPH_BACK_COLOR = clCream;
+  GRAPH_BACK_COLOR = $F0FBFF;  // clCream
-  GRAPH_WALL_COLOR = clGray;
+  GRAPH_WALL_COLOR = $808080;  // clGray
-  GRAPH_FLOOR_COLOR = clLtGray;
+  GRAPH_FLOOR_COLOR = $C0C0C0; // clLtGray
 
   TAB_FILE_FILTER = 'Tab field files (*.tab)|*.tab;*.TAB|Text files (*.txt)|*.txt;*.TXT|All files (*.*)|*.*';
   CSV_FILE_FILTER = 'Comma field files (*.csv)|*.csv;*.CSV|Text files (*.txt)|*.txt;*.TXT|All files (*.*)|*.*';

applications/lazstats/source/units/mathunit.pas

@@ -3,6 +3,7 @@ unit MathUnit;
 { extract some math functions from functionslib for easier testing }
 
 {$mode objfpc}{$H+}
+{$modeswitch advancedrecords}
 
 interface
 
@@ -20,15 +21,17 @@ function erfc(x: Double) : Double;
 
 function NormalDist(x: Double): Double;
 function NormalDistDensity(x, AMean, AStdDev: Double): Double;
+function InverseNormalDist(Probability: Double): Double;
 
 function Beta(a, b: Double): Extended;
 function BetaI(a,b,x: Double): Extended;
 
 function GammaLn(x: double): Extended;
 
-function tDist(x: Double; DF: Integer; OneSided: Boolean): Double;
+function tDist(x, DF: Double; OneSided: Boolean): Double;
-function tDensity(x: Double; DF: Integer): Double;
+function tDensity(x, DF: Double): Double;
-function ProbT(t, DF1: double): double;
+function ProbT(t, DF: double): double;
+function InverseT(Probability, DF: Double): Double;
 
 function FDensity(x: Double; DF1, DF2: Integer): Double;
 function ProbF(F, DF1, DF2: Double): Double;
@@ -44,7 +47,6 @@ function FactorialLn(n: Integer): Double;
 function PoissonPDF(n: integer; a: double): Double;
 function PoissonCDF(n: Integer; a: double): Double;
 
-
 procedure Calc_MaxMin(const AData: DblDyneVec; out AMax, AMin: Double);
 procedure Calc_MeanStdDev(const AData: DblDyneVec; out AMean, AStdDev: Double);
 procedure Calc_MeanVarStdDev(const AData: DblDyneVec; out AMean, AVariance, AStdDev: Double);
@@ -53,12 +55,38 @@ procedure Calc_SumSS(const AData: DblDyneVec; out Sum, SS: Double);
 
 function Calc_Median(const AData: DblDyneVec): Double;
 
+type
+  TBivariateRegressionResults = record
+    Slope, Intercept: Double;
+    XMean, YMean: Double;
+    XVariance, YVariance: Double;
+    XStdDev, YStdDev: Double;
+    StdErrorPredicted: Double;
+    R: Double;
+    t: Double;
+    SXX, SXY, SYY: Double;
+    Count, DF: Integer;
+    function ConfidenceLimits(x: Double; Upper: Boolean): Double;
+  end;
+
+procedure Calc_BivariateRegression(const xData, yData: DblDyneVec; AConfLevel: Double;
+  out AResults: TBivariateRegressionResults);
+
+procedure Exchange(var a, b: Double); overload;
+procedure Exchange(var a, b: Integer); overload;
+procedure Exchange(var a, b: String); overload;
+
+procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
+procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
+
+procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);   // not 100% tested...
+
 
 implementation
 
 uses
-  Math,
+  Math;
-  Utils;
+//  Utils;
 
 // Calculates the error function
 //                              /x
@@ -160,6 +188,80 @@ begin
 end;
 
 
+{ Obtains the inverse of the normal distribution (z), that is the argument for
+  the NormalDist() function to result in the given probability.
+
+  Probability = 0 ... 1 --> Result = -INF ... +INF
+
+  Algorithm by Peter John Acklam.
+  http://home.online.no/~pjacklam/notes/invnorm/index.html     }
+function InverseNormalDist(Probability: 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 Probability <= 0 then
+    Result := NegInfinity
+  else
+  if Probability < P_LOW then
+  begin
+    // rational approximation for lower region.
+    q := Sqrt(-2 * ln(Probability));
+    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 Probability <= P_HIGH then begin
+    // rational approximation for central region.
+    q := Probability - 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 Probability < 1 then begin
+    // rational approximation for upper region.
+    q := Sqrt(-2 * ln(1 - Probability));
+    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 := Infinity;
+end;
+
+
 function Beta(a, b: Double): Extended;
 begin
   if (a > 0) and (b > 0) then
@@ -268,14 +370,14 @@ begin
 end;
 
 // Calculates the (cumulative) t distribution function for DF degrees of freedom
-function tDist(x: Double; DF: Integer; OneSided: Boolean): Double;
+function tDist(x, DF: Double; OneSided: Boolean): Double;
 begin
   Result :=  betai(0.5*DF, 0.5, DF/(DF + sqr(x)));
   if OneSided then Result := Result * 0.5;
 end;
 
 // Returns the density curve for the t statistic with DF degrees of freedom
-function tDensity(x: Double; DF: Integer): Double;
+function tDensity(x, DF: Double): Double;
 var
   factor: Double;
 begin
@@ -285,12 +387,23 @@ end;
 
 { Returns the cumulative probability corresponding to a two-tailed t test with
   DF degrees of freedom. }
-function ProbT(t, DF1: double): double;
+function ProbT(t, DF: double): double;
 var
   F: double;
 begin
   F := t * t;
-  Result := ProbF(F, 1.0, DF1);
+  Result := ProbF(F, 1.0, DF);
+end;
+
+
+{ Returns the t value corresponding to a two-tailed t test probability. }
+function InverseT(Probability, DF: Double): double;
+var
+  z, w: double;
+begin
+  z := InverseNormalDist(Probability);
+  w := z * ((8.0 * DF + 3.0) / (1.0 + 8.0 * DF));
+  Result := sqrt(DF * (exp(w * w / DF) - 1.0));
 end;
 
 
@@ -590,6 +703,170 @@ begin
 end;
 
 
+// It is assumed that xData and yData contain at least 3 elements and
+// have the same count of elements.
+procedure Calc_BivariateRegression(const xData, yData: DblDyneVec; AConfLevel: Double;
+  out AResults: TBivariateRegressionResults);
+var
+  i: Integer;
+begin
+  with AResults do
+  begin
+    Count := Length(xData);
+
+    // Calculate means, variances, stddevs
+    Calc_MeanVarStdDevSS(xData, XMean, XVariance, XStdDev, SXX);
+    Calc_MeanVarStdDevSS(yData, YMean, YVariance, YStdDev, SYY);
+
+    SXY := 0;
+    for i := 0 to Count-1 do
+      SXY := SXY + xData[i] * yData[i];
+
+    R := (SXY - XMean * YMean * Count) / ((Count - 1) * XStdDev * YStdDev);
+    StdErrorPredicted := sqrt(1.0 - sqr(R)) * YStdDev * sqrt((Count - 1) / (Count - 2));
+
+    Slope := R * YStdDev / XStdDev;
+    Intercept := YMean - Slope * XMean;
+
+    DF := Count - 2;
+    t := InverseT(AConfLevel, DF);
+  end;
+end;
+
+
+function TBivariateRegressionResults.ConfidenceLimits(x: Double; Upper: Boolean): Double;
+var
+  yPred, seData: Double;
+begin
+  yPred := Intercept + Slope * x;
+  seData := StdErrorPredicted * sqrt(1.0 + 1/Count + sqr(x - XMean)/SXX);
+  if Upper then
+    Result := yPred + t*seData
+  else
+    Result := yPred - t*seData;
+end;
+
+
+procedure Exchange(var a, b: Double);
+var
+  tmp: Double;
+begin
+  tmp := a;
+  a := b;
+  b := tmp;
+end;
+
+procedure Exchange(var a, b: Integer);
+var
+  tmp: Integer;
+begin
+  tmp := a;
+  a := b;
+  b := tmp;
+end;
+
+procedure Exchange(var a, b: String);
+var
+  tmp: String;
+begin
+  tmp := a;
+  a := b;
+  b := tmp;
+end;
+
+procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
+var
+  i, j, N: Integer;
+begin
+  N := Length(X);
+  if (Y <> nil) and (N <> Length(Y)) then
+    raise Exception.Create('[SortOnX] Arrays must have the same length.');
+  if (Z <> nil) and (N <> Length(Z)) then
+    raise Exception.Create('[SortOnX] Arrays must have the same length.');
+
+  for i := 0 to N - 2 do
+  begin
+    for j := i + 1 to N - 1 do
+    begin
+      if X[i] > X[j] then //swap
+      begin
+        Exchange(X[i], X[j]);
+        if Y <> nil then
+          Exchange(Y[i], Y[j]);
+        if Z <> nil then
+          Exchange(Z[i], Z[j]);
+      end;
+    end;
+  end;
+end;
+
+// NOTE: The matrix Y is transposed relative to the typical usage in LazStats
+procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
+var
+  i, j, k, N, Ny: Integer;
+begin
+  N := Length(X);
+  if N <> Length(Y[0]) then
+    raise Exception.Create('[SortOnX] Arrays X and Y (2nd index) must have the same length');
+  Ny := Length(Y);
+
+  for i := 0 to N-2 do
+  begin
+    for j := i+1 to N-1 do
+      if X[i] > X[j] then
+      begin
+        Exchange(X[i], X[j]);
+        for k := 0 to Ny-1 do
+          Exchange(Y[k, i], Y[k, j]);
+      end;
+  end;
+end;
+
+procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
+
+  procedure DoQuickSort(L, R: Integer);
+  var
+    I,J: Integer;
+    P: Integer;
+  begin
+    repeat
+      I := L;
+      J := R;
+      P := (L+R) div 2;
+      repeat
+        while CompareValue(X[P], X[I]) > 0 do inc(I);
+        while CompareValue(X[P], X[J]) < 0 do dec(J);
+        if I <= J then begin
+          if I <> J then begin
+            Exchange(X[I], X[J]);
+            if Y <> nil then
+              Exchange(Y[I], Y[J]);
+            if Z <> nil then
+              Exchange(Z[I], Z[J]);
+          end;
+
+          if P = I then
+            P := J
+          else if P = J then
+            P := I;
+
+          inc(I);
+          dec(J);
+        end;
+      until I > J;
+
+      if L < J then
+        DoQuickSort(L, J);
+
+      L := I;
+    until I >= R;
+  end;
+
+begin
+  DoQuickSort(0, High(X));
+end;
+
+
 initialization
   InitFactLn();
 

applications/lazstats/source/units/utils.pas

@@ -21,15 +21,6 @@ function AnySelected(AListbox: TListBox): Boolean;
 procedure ErrorMsg(const AMsg: String);
 procedure ErrorMsg(const AMsg: String; const AParams: array of const);
 
-procedure Exchange(var a, b: Double); overload;
-procedure Exchange(var a, b: Integer); overload;
-procedure Exchange(var a, b: String); overload;
-
-procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
-procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
-
-procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);   // not 100% tested...
-
 function CenterString(S: String; Width: Integer): String;
 function IndexOfString(L: StrDyneVec; s: String): Integer;
 
@@ -107,125 +98,6 @@ begin
   ErrorMsg(Format(AMsg, AParams));
 end;
 
-procedure Exchange(var a, b: Double);
-var
-  tmp: Double;
-begin
-  tmp := a;
-  a := b;
-  b := tmp;
-end;
-
-procedure Exchange(var a, b: Integer);
-var
-  tmp: Integer;
-begin
-  tmp := a;
-  a := b;
-  b := tmp;
-end;
-
-procedure Exchange(var a, b: String);
-var
-  tmp: String;
-begin
-  tmp := a;
-  a := b;
-  b := tmp;
-end;
-
-procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
-var
-  i, j, N: Integer;
-begin
-  N := Length(X);
-  if (Y <> nil) and (N <> Length(Y)) then
-    raise Exception.Create('[SortOnX] Arrays must have the same length.');
-  if (Z <> nil) and (N <> Length(Z)) then
-    raise Exception.Create('[SortOnX] Arrays must have the same length.');
-
-  for i := 0 to N - 2 do
-  begin
-    for j := i + 1 to N - 1 do
-    begin
-      if X[i] > X[j] then //swap
-      begin
-        Exchange(X[i], X[j]);
-        if Y <> nil then
-          Exchange(Y[i], Y[j]);
-        if Z <> nil then
-          Exchange(Z[i], Z[j]);
-      end;
-    end;
-  end;
-end;
-
-// NOTE: The matrix Y is transposed relative to the typical usage in LazStats
-procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
-var
-  i, j, k, N, Ny: Integer;
-begin
-  N := Length(X);
-  if N <> Length(Y[0]) then
-    raise Exception.Create('[SortOnX] Arrays X and Y (2nd index) must have the same length');
-  Ny := Length(Y);
-
-  for i := 0 to N-2 do
-  begin
-    for j := i+1 to N-1 do
-      if X[i] > X[j] then
-      begin
-        Exchange(X[i], X[j]);
-        for k := 0 to Ny-1 do
-          Exchange(Y[k, i], Y[k, j]);
-      end;
-  end;
-end;
-
-procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
-
-  procedure DoQuickSort(L, R: Integer);
-  var
-    I,J: Integer;
-    P: Integer;
-  begin
-    repeat
-      I := L;
-      J := R;
-      P := (L+R) div 2;
-      repeat
-        while CompareValue(X[P], X[I]) > 0 do inc(I);
-        while CompareValue(X[P], X[J]) < 0 do dec(J);
-        if I <= J then begin
-          if I <> J then begin
-            Exchange(X[I], X[J]);
-            if Y <> nil then
-              Exchange(Y[I], Y[J]);
-            if Z <> nil then
-              Exchange(Z[I], Z[J]);
-          end;
-
-          if P = I then
-            P := J
-          else if P = J then
-            P := I;
-
-          inc(I);
-          dec(J);
-        end;
-      until I > J;
-
-      if L < J then
-        DoQuickSort(L, J);
-
-      L := I;
-    until I >= R;
-  end;
-
-begin
-  DoQuickSort(0, High(X));
-end;
-
 
 function CenterString(S: String; Width: Integer): String;
 var