From f386e5ab6bde35372d2a5b20d1659b9d49676612 Mon Sep 17 00:00:00 2001 From: wp Date: Tue, 21 Aug 2018 22:13:15 +0000 Subject: [PATCH] TAChart: Add support for confidence and prediction limits to TFitSeries. git-svn-id: trunk@58752 - --- components/tachart/demo/fit/Main.lfm | 90 ++++++++++++++++++++++------ components/tachart/demo/fit/Main.pas | 27 +++++++++ components/tachart/tafitlib.pas | 9 +++ components/tachart/tafitutils.pas | 55 ++++++++++++----- components/tachart/tafuncseries.pas | 56 +++++++++++++++-- 5 files changed, 199 insertions(+), 38 deletions(-) diff --git a/components/tachart/demo/fit/Main.lfm b/components/tachart/demo/fit/Main.lfm index ac48a66071..efde29ccc7 100644 --- a/components/tachart/demo/fit/Main.lfm +++ b/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 - Top = 452 + Height = 38 + 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 - Left = 132 + AnchorSideTop.Side = asrBottom + 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 diff --git a/components/tachart/demo/fit/Main.pas b/components/tachart/demo/fit/Main.pas index 5c52f2380c..07fe2322cc 100644 --- a/components/tachart/demo/fit/Main.pas +++ b/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.'); diff --git a/components/tachart/tafitlib.pas b/components/tachart/tafitlib.pas index b1ee7aad02..21843beb86 100644 --- a/components/tachart/tafitlib.pas +++ b/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); diff --git a/components/tachart/tafitutils.pas b/components/tachart/tafitutils.pas index b6d0cafc66..d10ff906a0 100644 --- a/components/tachart/tafitutils.pas +++ b/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; - aVarCovar: TArbFloatMatrix; ASignificanceLevel: Double = 0.05); + constructor Create(aFitResults: TFitResults; aAlpha: 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; - aVarCovar: TArbFloatMatrix; ASignificanceLevel: Double = 0.05); +constructor TFitStatistics.Create(aFitResults: TFitResults; + aAlpha: Double = 0.05); var i, j, L: Integer; begin - fN := aN; - fM := aM; - fSSR := aSSR; - fSSE := aSSE; - fSST := aSSR + aSSE; - fAlpha := ASignificanceLevel; - L := Length(aVarCovar); + fN := aFitResults.N; + fM := aFitResults.M; + fSSR := aFitResults.SSR; + fSSE := aFitResults.SSE; + fSST := aFitResults.SSR + aFitResults.SSE; + fAlpha := aAlpha; + 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. diff --git a/components/tachart/tafuncseries.pas b/components/tachart/tafuncseries.pas index db701243e8..d8f53cdf79 100644 --- a/components/tachart/tafuncseries.pas +++ b/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( - fitRes.N, fitRes.M, fitRes.SSR, fitRes.SSE, fitRes.CovarianceMatrix); + FFitStatistics := TFitStatistics.Create(fitRes, 1 - FConfidenceLevel); // 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;