TAChart: Add support for confidence and prediction limits to TFitSeries.

git-svn-id: trunk@58752 -

components/tachart/demo/fit/Main.lfm

@@ -552,7 +552,7 @@ object frmMain: TfrmMain
     TabOrder = 1
     object Chart: TChart
       Left = 0
-      Height = 452
+      Height = 449
       Top = 0
       Width = 560
       AxisList = <      
@@ -615,42 +615,73 @@ object frmMain: TfrmMain
         Source = ListChartSource
         OnFitComplete = FitCompleteHandler
       end
+      object UpperConfIntervalSeries: TFuncSeries
+        Active = False
+        Title = 'Confidence interval'
+        AxisIndexX = 1
+        AxisIndexY = 0
+        ExtentAutoY = True
+        Pen.Color = clRed
+        Pen.Style = psDot
+      end
+      object LowerConfIntervalSeries: TFuncSeries
+        Legend.Visible = False
+        Active = False
+        AxisIndexX = 1
+        AxisIndexY = 0
+        ExtentAutoY = True
+        Pen.Color = clRed
+        Pen.Style = psDot
+      end
+      object UpperPredIntervalSeries: TFuncSeries
+        Active = False
+        Title = 'Prediction interval'
+        AxisIndexX = 1
+        AxisIndexY = 0
+        ExtentAutoY = True
+        Pen.Color = clRed
+        Pen.Style = psDash
+      end
+      object LowerPredIntervalSeries: TFuncSeries
+        Legend.Visible = False
+        Active = False
+        AxisIndexX = 1
+        AxisIndexY = 0
+        ExtentAutoY = True
+        Pen.Color = clRed
+        Pen.Style = psDash
+      end
     end
     object pnlLog: TPanel
       Left = 0
-      Height = 35
+      Height = 38
-      Top = 452
+      Top = 449
       Width = 560
       Align = alBottom
       AutoSize = True
       BevelOuter = bvNone
-      ClientHeight = 35
+      ClientHeight = 38
       ClientWidth = 560
       TabOrder = 1
       object cbLogX: TCheckBox
         AnchorSideLeft.Control = pnlLog
         AnchorSideTop.Control = pnlLog
-        Left = 20
+        Left = 0
         Height = 19
-        Top = 8
+        Top = 0
         Width = 92
-        BorderSpacing.Left = 20
-        BorderSpacing.Top = 8
-        BorderSpacing.Bottom = 8
         Caption = 'Logarithmic x'
         OnClick = cbLogClick
         TabOrder = 0
       end
       object cbLogY: TCheckBox
         AnchorSideLeft.Control = cbLogX
-        AnchorSideLeft.Side = asrBottom
         AnchorSideTop.Control = cbLogX
-        AnchorSideTop.Side = asrCenter
+        AnchorSideTop.Side = asrBottom
-        Left = 132
+        Left = 0
         Height = 19
-        Top = 8
+        Top = 19
         Width = 93
-        BorderSpacing.Left = 20
         Caption = 'Logarithmic y'
         OnClick = cbLogClick
         TabOrder = 1
@@ -658,17 +689,42 @@ object frmMain: TfrmMain
       object cbShowErrorbars: TCheckBox
         AnchorSideLeft.Control = cbLogY
         AnchorSideLeft.Side = asrBottom
-        AnchorSideTop.Control = cbLogX
+        AnchorSideTop.Control = pnlLog
         AnchorSideTop.Side = asrCenter
-        Left = 245
+        Left = 113
         Height = 19
-        Top = 8
+        Top = 10
         Width = 102
         BorderSpacing.Left = 20
         Caption = 'Show error bars'
         OnChange = cbShowErrorbarsChange
         TabOrder = 2
       end
+      object cbShowConfidenceIntervals: TCheckBox
+        AnchorSideLeft.Control = cbShowErrorbars
+        AnchorSideLeft.Side = asrBottom
+        AnchorSideTop.Control = pnlLog
+        Left = 235
+        Height = 19
+        Top = 0
+        Width = 158
+        BorderSpacing.Left = 20
+        Caption = 'Show confidence intervals'
+        OnChange = cbShowConfidenceIntervalsChange
+        TabOrder = 3
+      end
+      object cbShowPredictionIntervals: TCheckBox
+        AnchorSideLeft.Control = cbShowConfidenceIntervals
+        AnchorSideTop.Control = cbShowConfidenceIntervals
+        AnchorSideTop.Side = asrBottom
+        Left = 235
+        Height = 19
+        Top = 19
+        Width = 153
+        Caption = 'Show prediction intervals'
+        OnChange = cbShowPredictionIntervalsChange
+        TabOrder = 4
+      end
     end
   end
   object Splitter1: TSplitter

components/tachart/demo/fit/Main.pas

@@ -22,6 +22,12 @@ type
     cbFitParam0Fixed: TCheckBox;
     cbFitParam1Fixed: TCheckBox;
     cbShowErrorbars: TCheckBox;
+    cbShowConfidenceIntervals: TCheckBox;
+    cbShowPredictionIntervals: TCheckBox;
+    UpperConfIntervalSeries: TFuncSeries;
+    LowerConfIntervalSeries: TFuncSeries;
+    UpperPredIntervalSeries: TFuncSeries;
+    LowerPredIntervalSeries: TFuncSeries;
     FitSeries: TFitSeries;
     cbFitRangeUseMin:TCheckBox;
     cbFitRangeUseMax:TCheckBox;
@@ -64,7 +70,9 @@ type
     procedure btnSaveClick(Sender: TObject);
     procedure cbDrawFitRangeOnlyClick(Sender: TObject);
     procedure cbFitEquationSelect(Sender: TObject);
+    procedure cbShowConfidenceIntervalsChange(Sender: TObject);
     procedure cbShowErrorbarsChange(Sender: TObject);
+    procedure cbShowPredictionIntervalsChange(Sender: TObject);
     procedure EdPointsCountChange(Sender: TObject);
     procedure FixedParamsChanged(Sender: TObject);
     procedure cbFitRangeUseMaxClick(Sender:TObject);
@@ -212,6 +220,12 @@ begin
   end;
 end;
 
+procedure TfrmMain.cbShowConfidenceIntervalsChange(Sender: TObject);
+begin
+  UpperConfIntervalSeries.Active := cbShowConfidenceIntervals.Checked;
+  LowerConfIntervalSeries.Active := cbShowConfidenceIntervals.Checked;
+end;
+
 procedure TfrmMain.cbShowErrorbarsChange(Sender: TObject);
 begin
   FitSeries.YErrorbars.Visible := cbShowErrorBars.Checked;
@@ -231,6 +245,12 @@ begin
     CreateData;
 end;
 
+procedure TfrmMain.cbShowPredictionIntervalsChange(Sender: TObject);
+begin
+  UpperPredIntervalSeries.Active := cbShowPredictionIntervals.Checked;
+  LowerPredIntervalSeries.Active := cbShowPredictionIntervals.Checked;
+end;
+
 procedure TfrmMain.EdPointsCountChange(Sender: TObject);
 begin
   CreateData;
@@ -448,6 +468,13 @@ begin
             Add('');
             Add('VARIANCE-COVARIANCE MATRIX');
             FitSeries.Statistics.Report_VarCovar(lbResults.Items);
+
+            {$IF FPC_FullVersion >= 30004}
+            UpperConfIntervalSeries.OnCalculate := @FitSeries.GetUpperConfidenceInterval;
+            LowerConfIntervalSeries.OnCalculate := @FitSeries.GetLowerConfidenceInterval;
+            UpperPredIntervalSeries.OnCalculate := @FitSeries.GetUpperPredictionInterval;
+            LowerPredIntervalSeries.OnCalculate := @FitSeries.GetLowerPredictionInterval;
+            {$IFEND}
           end;
         fitDimError:
           Add('The lengths of the data vectors do not match.');

components/tachart/tafitlib.pas

@@ -45,6 +45,8 @@ type
     M: Integer;                // Number of fit parameters
     SSR: ArbFloat;             // regression sum of squares (yhat - ybar)²
     SSE: ArbFloat;             // error sum of squares (yi - yhat)²
+    xbar: ArbFloat;            // mean x value
+    SSx: ArbFloat;             // sum of squares (xi - xbar)²
   end;
 
   { for compatibility with TAChart of Lazarus version <= 1.8.x }
@@ -351,6 +353,13 @@ begin
     end;
   end;
 
+  // Calculate x mean and sum of squares (needed for confidence intervals)
+  Result.xbar := 0;
+  for i := 0 to n - 1 do Result.xbar += x[i];
+  Result.xbar := Result.xbar / n;
+  Result.SSx := 0;
+  for i := 0 to n - 1 do Result.SSx += sqr(x[i] - Result.xbar);
+
   // Clean up
   SetLength(alpha, 0);
   SetLength(beta, 0);

components/tachart/tafitutils.pas

@@ -82,11 +82,15 @@ type
     fSSR: Double;       // regression sum of squares (yhat - ybar)^2
     fSSE: Double;       // error sum of squares (yi - yhat)^2
     fAlpha: Double;     // significance level for hypothesis tests
+    fxbar: Double;      // mean x value
+    fSSx: Double;       // sum of squares (xi - xbar)³
     fVarCovar: array of array of Double;
+    fTValue: Double;    // t-value
+    procedure CalcTValue;
     function GetVarCovar(i, j: Integer): Double;
+    procedure SetAlpha(AValue: Double);
   public
-    constructor Create(aN, aM: Integer; aSSR, aSSE: Double;
+    constructor Create(aFitResults: TFitResults; aAlpha: Double = 0.05);
-      aVarCovar: TArbFloatMatrix; ASignificanceLevel: Double = 0.05);
     procedure Report_ANOVA(AText: TStrings; ASeparator: String = ': ';
       ANumFormat: String = '%f');
     procedure Report_VarCovar(AText: TSTrings; ANumFormat: String = '%12.6f');
@@ -99,14 +103,19 @@ type
     property M: Integer read fM;
     function ReducedChi2: Double;
     function R2: Double;
-    function ResidualStandardError: Double;
+    function ResidualStdError: Double;
+    property Alpha: Double read FAlpha write SetAlpha;
     property SST: Double read fSST;
     property SSR: Double read fSSR;
     property SSE: Double read fSSE;
     property VarCovar[i, j: Integer]: Double read GetVarCovar;
+    property xBar: Double read fXBar;
+    property SSx: Double read fSSx;
+  public
     {$IF FPC_FullVersion >= 30004}
     function Fcrit: Double;
     function pValue: Double;
+    property tValue: Double read ftValue;
     {$ENDIF}
   end;
 
@@ -314,22 +323,25 @@ end;
 
 { TFitStatistics }
 
-constructor TFitStatistics.Create(aN, aM: Integer; aSSR, aSSE: Double;
+constructor TFitStatistics.Create(aFitResults: TFitResults;
-  aVarCovar: TArbFloatMatrix; ASignificanceLevel: Double = 0.05);
+  aAlpha: Double = 0.05);
 var
   i, j, L: Integer;
 begin
-  fN := aN;
+  fN := aFitResults.N;
-  fM := aM;
+  fM := aFitResults.M;
-  fSSR := aSSR;
+  fSSR := aFitResults.SSR;
-  fSSE := aSSE;
+  fSSE := aFitResults.SSE;
-  fSST := aSSR + aSSE;
+  fSST := aFitResults.SSR + aFitResults.SSE;
-  fAlpha := ASignificanceLevel;
+  fAlpha := aAlpha;
-  L := Length(aVarCovar);
+  L := Length(aFitResults.CovarianceMatrix);
   SetLength(fVarCovar, L, L);
   for j := 0 to L-1 do
     for i := 0 to L-1 do
-      fVarCovar[i, j] := aVarCovar[i, j];
+      fVarCovar[i, j] := aFitResults.CovarianceMatrix[i, j];
+  fXBar := aFitResults.XBar;
+  fSSx := aFitResults.SSx;
+  CalcTValue;
 end;
 
 { Coefficient of determination, adjusted to number of data points and fit
@@ -339,6 +351,13 @@ begin
   Result := 1.0 - (1.0 - R2) * (N - 1) / DOF;
 end;
 
+procedure TFitStatistics.CalcTValue;
+begin
+  {$IF FPC_FullVersion >= 30004}
+  fTValue := invtdist(fAlpha, fN - fM, 2)
+  {$IFEND}
+end;
+
 { Total variance of data values minus calculated values, weighted by
   data error.
   For a "moderately" good fit Chi2 is approximately equal to the degrees of
@@ -412,7 +431,7 @@ begin
 end;
 
 { Mean residual standard error of fit: The smaller the better }
-function TFitStatistics.ResidualStandardError: Double;
+function TFitStatistics.ResidualStdError: Double;
 begin
   if DOF > 0 then
     Result := sqrt(SSE / DOF)
@@ -438,7 +457,7 @@ begin
   AText.Add(rsFitAdjCoefficientOfDetermination + ASeparator + Format(ANumFormat, [AdjR2]));
   AText.Add(rsFitChiSquared + ASeparator + Format(ANumFormat, [Chi2]));
   AText.Add(rsFitReducedChiSquared + ASeparator + Format(ANumFormat, [ReducedChi2]));
-  AText.Add(rsFitResidualStandardError + ASeparator + Format(ANumFormat, [ResidualStandardError]));
+  AText.Add(rsFitResidualStandardError + ASeparator + Format(ANumFormat, [ResidualStdError]));
   AText.Add(rsFitVarianceRatio + ASeparator + Format(ANumFormat, [F]));
   {
   AText.Add(Format('Fcrit(%d, %d)', [M-1, DOF]) + ASeparator +
@@ -471,5 +490,11 @@ begin
   end;
 end;
 
+procedure TFitStatistics.SetAlpha(AValue: Double);
+begin
+  fAlpha := AValue;
+  CalcTValue;
+end;
+
 end.
 

components/tachart/tafuncseries.pas

@@ -301,9 +301,6 @@ type
     function GetParam(AIndex: Integer): Double;
     function GetParamCount: Integer;
     function GetParamError(AIndex: Integer): Double;
-    {$IF FPC_FullVersion >= 30004}
-    function GetParam_pValue(AIndex: Integer): Double;
-    {$IFEND}
     function GetParam_RawError(AIndex: Integer): Double;
     function GetParam_RawValue(AIndex: Integer): Double;
     function GetParam_tValue(AIndex: Integer): Double;
@@ -316,6 +313,10 @@ type
     procedure SetParamCount(AValue: Integer);
     procedure SetPen(AValue: TChartPen);
     procedure SetStep(AValue: TFuncSeriesStep);
+    {$IF FPC_FullVersion >= 30004}
+    procedure GetInterval(const Ax: Double; out AY: Double; IsUpper, IsPrediction: Boolean);
+    function GetParam_pValue(AIndex: Integer): Double;
+    {$IFEND}
   strict protected
     procedure CalcXRange(out AXMin, AXMax: Double);
     function TransformX(AX: Double): Extended; inline;
@@ -337,6 +338,10 @@ type
     function FitParams: TDoubleDynArray;
     {$IF FPC_FullVersion >= 30004}
     procedure GetConfidenceLimits(AIndex: Integer; out ALower, AUpper: Double);
+    procedure GetLowerConfidenceInterval(const Ax: Double; out AY: Double);
+    procedure GetUpperConfidenceInterval(const Ax: Double; out AY: Double);
+    procedure GetLowerPredictionInterval(const Ax: Double; out AY: Double);
+    procedure GetUpperPredictionInterval(const Ax: Double; out AY: Double);
     {$IFEND}
     function GetFitEquationString(
       ANumFormat: String; AXText: String = 'x'; AYText: String = 'y'): String;
@@ -1696,8 +1701,7 @@ var
 
     // Analysis of variance, variance-covariance matrix
     FFitStatistics.Free;
-    FFitStatistics := TFitStatistics.Create(
+    FFitStatistics := TFitStatistics.Create(fitRes, 1 - FConfidenceLevel);
-      fitRes.N, fitRes.M, fitRes.SSR, fitRes.SSE, fitRes.CovarianceMatrix);
 
     // State of the fit
     FState := fpsValid;
@@ -1733,7 +1737,8 @@ begin
   val := GetParam_RawValue(AIndex);
   sig := GetParam_RawError(AIndex);
   alpha := 1.0 - FConfidenceLevel;
-  t := invtdist(alpha, Statistics.DOF, 2);
+  t := Statistics.tValue;
+//  t := invtdist(alpha, Statistics.DOF, 2);
   ALower := val - sig*t;
   AUpper := val + sig*t;
   if (FFitEquation in [feExp, fePower]) and (AIndex = 0) then begin
@@ -1741,6 +1746,45 @@ begin
     AUpper := exp(AUpper);
   end;
 end;
+
+procedure TFitSeries.GetInterval(const aX: Double; out AY: Double;
+  IsUpper, IsPrediction: Boolean);
+var
+  y: Double;
+  dy: Double;
+  Offs: Double;
+begin
+  offs := IfThen(IsPrediction, 1, 0);
+  with Statistics do begin
+    y := TransformY(Calculate(AX));
+    dy := tValue * ResidualStdError * sqrt(offs + 1/N + sqr(AX - xBar) / SSx);
+    if IsUpper then
+      AY := y + dy
+    else
+      AY := y - dy;
+    if (FFitEquation in [feExp, fePower]) then AY := exp(AY);
+  end;
+end;
+
+procedure TFitSeries.GetLowerConfidenceInterval(const AX: Double; out AY: Double);
+begin
+  GetInterval(AX, AY, false, false);
+end;
+
+procedure TFitSeries.GetUpperConfidenceInterval(const AX: Double; out AY: Double);
+begin
+  GetInterval(AX, AY, true, false);
+end;
+
+procedure TFitSeries.GetLowerPredictionInterval(const AX: Double; out AY: Double);
+begin
+  GetInterval(AX, AY, false, true);
+end;
+
+procedure TFitSeries.GetUpperPredictionInterval(const AX: Double; out AY: Double);
+begin
+  GetInterval(AX, AY, true, true);
+end;
 {$IFEND}
 
 function TFitSeries.GetFitEquationString(ANumFormat: String; AXText: String;