1 {
2
3 *****************************************************************************
4 See the file COPYING.modifiedLGPL.txt, included in this distribution,
5 for details about the license.
6 *****************************************************************************
7
8 Authors: Alexander Klenin
9
10 }
11 unit TAFitUtils;
12
13 {$H+}
14
15 interface
16
17 uses
18 SysUtils, Classes, TAChartUtils, TAFitLib;
19
20 type
21 TFitEquation = (
22 fePolynomial, // y = b0 + b1*x + b2*x^2 + ... bn*x^n
23 feLinear, // y = a + b*x
24 feExp, // y = a * exp(b * x)
25 fePower, // y = a * x^b
26 feCustom // y = b0 + b1*F1(x) + b2*F2(x) + ... bn*Fn(x),
providednull27 // Fi(x) = custom "fit base function" provided by calling SetFitBasisFunc() method
28 );
29
30 IFitEquationText = interface
BasisFuncsnull31 function BasisFuncs(const ATexts: array of string): IFitEquationText;
DecimalSeparatornull32 function DecimalSeparator(AValue: Char): IFitEquationText;
Equationnull33 function Equation(AEquation: TFitEquation): IFitEquationText;
Xnull34 function X(AText: String): IFitEquationText;
Ynull35 function Y(AText: String): IFitEquationText;
NumFormatnull36 function NumFormat(AFormat: String): IFitEquationText;
NumFormatsnull37 function NumFormats(const AFormats: array of String): IFitEquationText;
Paramsnull38 function Params(const AParams: array of Double): IFitEquationText;
TextFormatnull39 function TextFormat(AFormat: TChartTextFormat): IFitEquationText;
Getnull40 function Get: String;
41 end;
42
43 TFitEmptyEquationText = class(TInterfacedObject, IFitEquationText)
44 public
BasisFuncsnull45 function BasisFuncs(const ATexts: array of string): IFitEquationText;
DecimalSeparatornull46 function DecimalSeparator(AValue: Char): IFitEquationText;
Equationnull47 function Equation(AEquation: TFitEquation): IFitEquationText;
Xnull48 function X(AText: String): IFitEquationText;
Ynull49 function Y(AText: String): IFitEquationText;
NumFormatnull50 function NumFormat(AFormat: String): IFitEquationText;
NumFormatsnull51 function NumFormats(const AFormats: array of String): IFitEquationText;
Paramsnull52 function Params(const AParams: array of Double): IFitEquationText;
TextFormatnull53 function TextFormat(AFormat: TChartTextFormat): IFitEquationText;
Getnull54 function Get: String;
55 end;
56
57 TFitEquationText = class(TInterfacedObject, IFitEquationText)
58 strict private
59 FBasisFunc: array of string;
60 FDecSep: Char;
61 FEquation: TFitEquation;
62 FX: String;
63 FY: String;
64 FNumFormat: String;
65 FNumFormats: array of String;
66 FParams: array of Double;
67 FTextFormat: TChartTextFormat;
GetNumFormatnull68 function GetNumFormat(AIndex: Integer): String;
69 public
70 constructor Create;
BasisFuncsnull71 function BasisFuncs(const ATexts: array of string): IFitEquationText;
DecimalSeparatornull72 function DecimalSeparator(AValue: Char): IFitEquationText;
Equationnull73 function Equation(AEquation: TFitEquation): IFitEquationText;
Xnull74 function X(AText: String): IFitEquationText;
Ynull75 function Y(AText: String): IFitEquationText;
NumFormatnull76 function NumFormat(AFormat: String): IFitEquationText;
NumFormatsnull77 function NumFormats(const AFormats: array of String): IFitEquationText;
Paramsnull78 function Params(const AParams: array of Double): IFitEquationText;
TextFormatnull79 function TextFormat(AFormat: TChartTextFormat): IFitEquationText;
Getnull80 function Get: String;
81 end;
82
83 TFitStatistics = class(TObject)
84 private
85 fN: Integer; // number of observations
86 fM: Integer; // number of (adjusted) fit parameters (fixed params not included)
87 fSST: Double; // total sum of squares (yi - ybar)^2
88 fSSR: Double; // regression sum of squares (yhat - ybar)^2
89 fSSE: Double; // error sum of squares (yi - yhat)^2
90 fAlpha: Double; // significance level for hypothesis tests
91 fxbar: Double; // mean x value
92 fSSx: Double; // sum of squares (xi - xbar)³
93 fVarCovar: array of array of Double;
94 fTValue: Double; // t-value
95 procedure CalcTValue;
GetVarCovarnull96 function GetVarCovar(i, j: Integer): Double;
97 procedure SetAlpha(AValue: Double);
98 public
99 constructor Create(aFitResults: TFitResults; aAlpha: Double = 0.05);
100 procedure Report_ANOVA(AText: TStrings; ASeparator: String = ': ';
101 ANumFormat: String = '%f'; AExpFormat: String = '%.3e'; NaNStr: String = 'n/a');
102 procedure Report_VarCovar(AText: TSTrings; ANumFormat: String = '%12.6f');
103 public
AdjR2null104 function AdjR2: Double;
Chi2null105 function Chi2: Double;
DOFnull106 function DOF: Integer; // Degrees of freedom
Fnull107 function F: Double;
108 property N: Integer read fN; // Number of data points
109 property M: Integer read fM; // Number of fit parameters
ReducedChi2null110 function ReducedChi2: Double;
R2null111 function R2: Double;
ResidualStdErrornull112 function ResidualStdError: Double;
113 property Alpha: Double read FAlpha write SetAlpha;
114 property SST: Double read fSST;
115 property SSR: Double read fSSR;
116 property SSE: Double read fSSE;
117 property VarCovar[i, j: Integer]: Double read GetVarCovar;
118 property xBar: Double read fXBar;
119 property SSx: Double read fSSx;
120 public
121 {$IF FPC_FullVersion >= 30004}
Fcritnull122 function Fcrit: Double;
pValuenull123 function pValue: Double;
124 property tValue: Double read ftValue;
125 {$ENDIF}
126 end;
127
128 operator := (AEq: IFitEquationText): String; inline;
129
130 implementation
131
132 uses
133 Math, spe, TAChartStrConsts;
134
135 operator := (AEq: IFitEquationText): String;
136 begin
137 Result := AEq.Get;
138 end;
139
140
141 { TFitEmptyEquationText }
142
BasisFuncsnull143 function TFitEmptyEquationText.BasisFuncs(const ATexts: array of string): IFitEquationText;
144 begin
145 Unused(ATexts);
146 Result := Self;
147 end;
148
DecimalSeparatornull149 function TFitEmptyEquationText.DecimalSeparator(AValue: Char): IFitEquationText;
150 begin
151 Unused(AValue);
152 Result := Self;
153 end;
154
TFitEmptyEquationText.Equationnull155 function TFitEmptyEquationText.Equation(
156 AEquation: TFitEquation): IFitEquationText;
157 begin
158 Unused(AEquation);
159 Result := Self;
160 end;
161
TFitEmptyEquationText.Getnull162 function TFitEmptyEquationText.Get: String;
163 begin
164 Result := '';
165 end;
166
NumFormatnull167 function TFitEmptyEquationText.NumFormat(AFormat: String): IFitEquationText;
168 begin
169 Unused(AFormat);
170 Result := Self;
171 end;
172
TFitEmptyEquationText.NumFormatsnull173 function TFitEmptyEquationText.NumFormats(
174 const AFormats: array of String): IFitEquationText;
175 begin
176 Unused(AFormats);
177 Result := Self;
178 end;
179
Paramsnull180 function TFitEmptyEquationText.Params(
181 const AParams: array of Double): IFitEquationText;
182 begin
183 Unused(AParams);
184 Result := Self;
185 end;
186
TFitEmptyEquationText.TextFormatnull187 function TFitEmptyEquationText.TextFormat(AFormat: TChartTextFormat): IFitEquationText;
188 begin
189 Unused(AFormat);
190 Result := Self;
191 end;
192
Xnull193 function TFitEmptyEquationText.X(AText: String): IFitEquationText;
194 begin
195 Unused(AText);
196 Result := Self;
197 end;
198
Ynull199 function TFitEmptyEquationText.Y(AText: String): IFitEquationText;
200 begin
201 Unused(AText);
202 Result := Self;
203 end;
204
205 { TFitEquationText }
206
207 constructor TFitEquationText.Create;
208 begin
209 FX := 'x';
210 FY := 'y';
211 FNumFormat := '%.9g';
212 FDecSep := DefaultFormatSettings.DecimalSeparator;
213 end;
214
BasisFuncsnull215 function TFitEquationText.BasisFuncs(const ATexts: array of string): IFitEquationText;
216 var
217 i: Integer;
218 begin
219 SetLength(FBasisFunc, Length(ATexts));
220 for i := 0 to High(FBasisFunc) do
221 FBasisFunc[i] := ATexts[i];
222 Result := Self;
223 end;
224
DecimalSeparatornull225 function TFitEquationText.DecimalSeparator(AValue: Char): IFitEquationText;
226 begin
227 FDecSep := AValue;
228 Result := self;
229 end;
230
TFitEquationText.Equationnull231 function TFitEquationText.Equation(AEquation: TFitEquation): IFitEquationText;
232 begin
233 FEquation := AEquation;
234 Result := Self;
235 end;
236
TFitEquationText.Getnull237 function TFitEquationText.Get: String;
238 const
239 TIMES: array[boolean] of string = ('*', '·');
240 var
241 fs: TFormatSettings;
242 res: String;
243
termnull244 // Returns the function term, e.g. "x^3"
245 function FuncTerm(i: Integer): String;
246 const
247 POWER: array[boolean] of string = ('%s^%d', '%s<sup>%d</sup>');
248 begin
249 if FEquation in [feCustom] then
250 Result := FBasisFunc[i]
251 else
252 if i = 0 then
253 Result := ''
254 else
255 if i = 1 then
256 Result := FX
257 else
258 Result := Format(POWER[FTextFormat = tfHTML], [FX, i]);
259 end;
260
261 // Creates a product term "value*f(x)".
262 // "value*" is omitted if 1. "*f(x)" is omitted if constant.
ProductTermnull263 function ProductTerm(i: Integer): String;
264 var
265 fx: String;
266 begin
267 if FParams[i] = 0.0 then
268 exit('');
269 fx := FuncTerm(i);
270 if abs(FParams[i]) <> 1.0 then begin
271 Result := Format(GetNumFormat(i), [FParams[i]], fs);
272 if fx <> '' then
273 Result := Result + TIMES[FTextFormat = tfHTML] + fx;
274 end else
275 if FParams[i] = -1.0 then
276 Result := '-' + FX
277 else
278 Result := FX;
279 end;
280
281 // First term in expression: no plus sign
282 // Other terms: both plus and minus signs, enclosed by spaces.
FixSignnull283 function FixSign(s: String): String;
284 begin
285 if (s <> '') and (res <> '') then begin
286 if s[1] = '-' then begin
287 Insert(' ', s, 2);
288 s := ' ' + s;
289 end else
290 s := ' + ' + s;
291 end;
292 Result := s;
293 end;
294
295 // Creates the constant term for feExp or fePower; includes the multiplication sign.
ConstTermnull296 function ConstTerm: String;
297 begin
298 if FParams[0] = 1.0 then
299 Result := ''
300 else if FParams[0] = -1.0 then
301 Result := '-'
302 else
303 Result := Format(GetNumFormat(0), [FParams[0]], fs) + TIMES[FTextFormat = tfHTML];
304 end;
305
306 // Creates the exponential term, e.g. 'exp(-1.2*x)' or 'e<sup>-1.2·x</sup>'
ExpTermnull307 function ExpTerm: string;
308 const
309 EXP: array[boolean] of String = ('exp(%s)', 'e<sup>%s</sup>');
310 begin
311 Result := Format(EXP[FTextFormat = tfHTML], [ProductTerm(1)]);
312 end;
313
314 // Creates the power term, e.g. 'x^1.2' or 'x<sup>1.2</sup>'
PowerTermnull315 function PowerTerm: String;
316 var
317 mask: String;
318 begin
319 if FTextFormat = tfNormal then
320 mask := '%s^' + GetNumFormat(1)
321 else
322 mask := '%s<sup>' + GetNumFormat(1) + '</sup>';
323 Result := Format(mask, [FX, FParams[1]], fs);
324 end;
325
326 var
327 i: Integer;
328 begin
329 if Length(FParams) = 0 then
330 exit('');
331
332 fs := DefaultFormatSettings;
333 fs.DecimalSeparator := FDecSep;
334
335 Result := Format('%s = ', [FY]);
336 case FEquation of
337 feLinear, fePolynomial, feCustom:
338 begin
339 res := '';
340 for i := 0 to High(FParams) do
341 res += FixSign(ProductTerm(i));
342 Result += res;
343 end;
344 feExp:
345 Result += ConstTerm + ExpTerm;
346 fePower:
347 Result += ConstTerm + PowerTerm;
348 end;
349 end;
350
GetNumFormatnull351 function TFitEquationText.GetNumFormat(AIndex: Integer): String;
352 begin
353 if AIndex < Length(FNumFormats) then
354 Result := FNumFormats[AIndex]
355 else
356 Result := FNumFormat;
357 end;
358
TFitEquationText.NumFormatnull359 function TFitEquationText.NumFormat(AFormat: String): IFitEquationText;
360 begin
361 FNumFormat := AFormat;
362 Result := Self;
363 end;
364
TFitEquationText.NumFormatsnull365 function TFitEquationText.NumFormats(
366 const AFormats: array of String): IFitEquationText;
367 var
368 i: Integer;
369 begin
370 SetLength(FNumFormats, Length(AFormats));
371 for i := 0 to High(AFormats) do
372 FNumFormats[i] := AFormats[i];
373 Result := Self;
374 end;
375
TFitEquationText.Paramsnull376 function TFitEquationText.Params(
377 const AParams: array of Double): IFitEquationText;
378 var
379 i: Integer;
380 begin
381 SetLength(FParams, Length(AParams));
382 for i := 0 to High(AParams) do
383 FParams[i] := AParams[i];
384 Result := Self;
385 end;
386
TFitEquationText.TextFormatnull387 function TFitEquationText.TextFormat(AFormat: TChartTextFormat): IFitEquationText;
388 begin
389 FTextFormat := AFormat;
390 Result := Self;
391 end;
392
Xnull393 function TFitEquationText.X(AText: String): IFitEquationText;
394 begin
395 FX := AText;
396 Result := Self;
397 end;
398
Ynull399 function TFitEquationText.Y(AText: String): IFitEquationText;
400 begin
401 FY := AText;
402 Result := Self;
403 end;
404
405
406 { TFitStatistics }
407
408 constructor TFitStatistics.Create(aFitResults: TFitResults;
409 aAlpha: Double = 0.05);
410 var
411 i, j, L: Integer;
412 begin
413 fN := aFitResults.N;
414 fM := aFitResults.M;
415 fSSR := aFitResults.SSR;
416 fSSE := aFitResults.SSE;
417 fSST := aFitResults.SSR + aFitResults.SSE;
418 fAlpha := aAlpha;
419 L := Length(aFitResults.CovarianceMatrix);
420 SetLength(fVarCovar, L, L);
421 for j := 0 to L-1 do
422 for i := 0 to L-1 do
423 fVarCovar[i, j] := aFitResults.CovarianceMatrix[i, j];
424 fXBar := aFitResults.XBar;
425 fSSx := aFitResults.SSx;
426 CalcTValue;
427 end;
428
429 { Coefficient of determination, adjusted to number of data points and fit
430 parameters. Should be close to 1 ("good" fit). "0" means: "poor" fit }
AdjR2null431 function TFitStatistics.AdjR2: Double;
432 begin
433 if DOF > 0 then
434 Result := 1.0 - (1.0 - R2) * (N - 1) / DOF
435 else
436 Result := NaN;
437 end;
438
439 procedure TFitStatistics.CalcTValue;
440 begin
441 fTValue := NaN;
442 {$IF FPC_FullVersion >= 30004}
443 if (fAlpha > 0) and (fN > fM) then
444 fTValue := invtdist(fAlpha, fN - fM, 2)
445 {$IFEND}
446 end;
447
448 { Total variance of data values minus calculated values, weighted by
449 data error.
450 For a "moderately" good fit Chi2 is approximately equal to the degrees of
451 freedom (DOF). }
TFitStatistics.Chi2null452 function TFitStatistics.Chi2: Double;
453 begin
454 Result := SSE;
455 end;
456
457 { Degrees of freedom }
DOFnull458 function TFitStatistics.DOF: Integer;
459 begin
460 Result := N - M;
461 end;
462
Fnull463 function TFitStatistics.F: Double;
464 begin
465 if (M > 1) and (N <> M) and (SSE <> 0) then
466 Result := (SSR/(M-1)) / (SSE/(N-M))
467 else
468 Result := NaN;
469 end;
470
471 {$IF FPC_FullVersion >= 30004}
Fcritnull472 function TFitStatistics.Fcrit: Double;
473 begin
474 if (M = 1) then
475 Result := NaN
476 else
477 Result := InvFDist(FAlpha, M-1, N-M);
478 end;
479 {$IFEND}
480
GetVarCovarnull481 function TFitStatistics.GetVarCovar(i, j: Integer): Double;
482 begin
483 Result := fVarCovar[i, j];
484 end;
485
486 {$IF FPC_FullVersion >= 30004}
487 { Probability that the scatter of the data around the fitted curve is by chance.
488 Should be several 0.1, the higher the better.
489 According to Numerical Recipes, very small (<< 0.1) values mean
490 - wrong model
491 - error bars too small
492 - errors are not normally distributed
493 Values close to 1.0 mean
494 - error bars overestimaged. }
TFitStatistics.pValuenull495 function TFitStatistics.pValue: Double;
496 begin
497 if DOF > 0 then
498 Result := chi2dist(Chi2, DOF)
499 else
500 Result := NaN;
501 end;
502 {$IFEND}
503
504 { Variance normalized to the degrees of freedem. Should be about 1 for
505 a "moderately" good fit. }
TFitStatistics.ReducedChi2null506 function TFitStatistics.ReducedChi2: Double;
507 begin
508 if (DOF > 0) then
509 Result := SSE / DOF
510 else
511 Result := NaN;
512 end;
513
514 { Coefficient of determination: 0 -> "poor" fit, 1 -> "good" fit }
R2null515 function TFitStatistics.R2: Double;
516 begin
517 Result := SSR / SST;
518 end;
519
520 { Mean residual standard error of fit: The smaller the better }
ResidualStdErrornull521 function TFitStatistics.ResidualStdError: Double;
522 begin
523 if DOF > 0 then
524 Result := sqrt(SSE / DOF)
525 else
526 Result := NaN;
527 end;
528
529 procedure TFitStatistics.Report_ANOVA(AText: TStrings; ASeparator: String = ': ';
530 ANumFormat: String = '%f'; AExpFormat: String = '%.3e'; NaNStr: String = 'n/a');
531 const
532 PRECISION = 3;
533 begin
534 AText.Add(rsFitNumObservations + ASeparator + IntToStr(N));
535 AText.Add(rsFitNumFitParams + ASeparator + IntToStr(M));
536 AText.Add(rsFitDegreesOfFreedom + ASeparator + IntToStr(DOF));
537 AText.Add(rsFitTotalSumOfSquares + ASeparator + FloatToStrEx(SST, PRECISION, ANumFormat, AExpFormat, NaNStr));
538 AText.Add(rsFitRegressionSumOfSquares + ASeparator + FloatToStrEx(SSR, PRECISION, ANumFormat, AExpFormat, NaNStr));
539 AText.Add(rsFitErrorSumOfSquares + ASeparator + FloatToStrEx(SSE, PRECISION, ANumFormat, AExpFormat, NaNStr));
540 AText.Add(rsFitCoefficientOfDetermination + ASeparator + FloatToStrEx(R2, PRECISION, ANumFormat, '', NaNStr));
541 AText.Add(rsFitAdjCoefficientOfDetermination + ASeparator + FloatToStrEx(AdjR2, PRECISION, ANumFormat, '', NaNStr));
542 AText.Add(rsFitChiSquared + ASeparator + FloatToStrEx(Chi2, PRECISION, ANumFormat, AExpFormat, NaNStr));
543 AText.Add(rsFitReducedChiSquared + ASeparator + FloatToStrEx(ReducedChi2, PRECISION, ANumFormat, AExpFormat, NaNStr));
544 AText.Add(rsFitResidualStandardError + ASeparator + FloatToStrEx(ResidualStdError, PRECISION, ANumFormat, AExpFormat, NaNStr));
545 AText.Add(rsFitVarianceRatio + ASeparator + FloatToStrEx(F, PRECISION, ANumFormat, AExpFormat, NaNStr));
546 {
547 AText.Add(Format('Fcrit(%d, %d)', [M-1, DOF]) + ASeparator +
548 Format(IfThen(Fcrit < 1E-3, FMT, ANumFormat), [Fcrit]));
549 }
550 {$IF FPC_FullVersion >= 30004}
551 AText.Add(rsFitTValue + ASeparator + FloatToStrEx(FtValue, PRECISION, ANumFormat, AExpFormat, NaNStr));
552 AText.Add(rsFitPValue + ASeparator + FloatToStrEx(pValue, PRECISION, ANumFormat, AExpFormat, NaNStr));
553 {$IFEND}
554 end;
555
556 procedure TFitStatistics.Report_VarCovar(AText: TStrings; ANumFormat: String = '%12.6f');
557 var
558 i, j: Integer;
559 s, t: String;
560 w: Integer;
561 begin
562 t := Copy(ANumFormat, 1, pos('.', ANumFormat)-1);
563 Delete(t, 1, 1);
564 w := StrToInt(t);
565
566 for i := 0 to M-1 do begin
567 s := '';
568 for j := 0 to M-1 do
569 if IsNaN(VarCovar[i, j]) then
570 s := s + Format('%*s', [w, '---'])
571 else
572 s := s + Format(ANumFormat, [VarCovar[i, j]]);
573 AText.Add(s);
574 end;
575 end;
576
577 procedure TFitStatistics.SetAlpha(AValue: Double);
578 begin
579 fAlpha := AValue;
580 CalcTValue;
581 end;
582
583 end.
584
585