diff --git a/applications/lazstats/source/LazStats.lpi b/applications/lazstats/source/LazStats.lpi index c5e6015c7..07a93c1d3 100644 --- a/applications/lazstats/source/LazStats.lpi +++ b/applications/lazstats/source/LazStats.lpi @@ -49,7 +49,7 @@ - + @@ -1400,6 +1400,11 @@ + + + + + diff --git a/applications/lazstats/source/LazStats.lpr b/applications/lazstats/source/LazStats.lpr index 20ca88a54..a284fbee6 100644 --- a/applications/lazstats/source/LazStats.lpr +++ b/applications/lazstats/source/LazStats.lpr @@ -7,8 +7,8 @@ uses cthreads, {$ENDIF}{$ENDIF} Interfaces, // this includes the LCL widgetset - Forms, tachartlazaruspkg, tachartprint, lhelpcontrolpkg, - Globals, LicenseUnit, OptionsUnit, MainDM, MainUnit; //, utils, chartunit; + Forms, tachartlazaruspkg, tachartprint, lhelpcontrolpkg, Globals, LicenseUnit, + OptionsUnit, MainDM, MainUnit, MathUnit; //, utils, chartunit; {$R LazStats.res} diff --git a/applications/lazstats/source/forms/analysis/statistical_process_control/sigmachartunit.pas b/applications/lazstats/source/forms/analysis/statistical_process_control/sigmachartunit.pas index 18f6cedd4..4048c264c 100644 --- a/applications/lazstats/source/forms/analysis/statistical_process_control/sigmachartunit.pas +++ b/applications/lazstats/source/forms/analysis/statistical_process_control/sigmachartunit.pas @@ -51,7 +51,8 @@ var implementation uses - Math; + Math, + MathUnit; { TSigmaChartFrm } @@ -252,9 +253,9 @@ begin D3Value := D3[grpsize-1]; D4Value := D4[grpsize-1]; C4 := sqrt(2.0 / (grpsize - 1)); - gamma := exp(gammln(grpsize / 2.0)); + gamma := exp(GammaLn(grpsize / 2.0)); C4 := C4 * gamma; - gamma := exp(gammln((grpsize-1) / 2.0)); + gamma := exp(GammaLn((grpsize-1) / 2.0)); C4 := C4 / gamma; B := GrandSigma * sqrt(1.0 - (C4 * C4)) / C4; UCL := GrandSigma + 3.0 * B; diff --git a/applications/lazstats/source/forms/simulations/distribunit.lfm b/applications/lazstats/source/forms/simulations/distribunit.lfm index c83831814..7e49c860e 100644 --- a/applications/lazstats/source/forms/simulations/distribunit.lfm +++ b/applications/lazstats/source/forms/simulations/distribunit.lfm @@ -122,7 +122,7 @@ object DistribFrm: TDistribFrm Minors = <> Title.LabelFont.Style = [fsBold] Title.Visible = True - Title.Caption = 'Scale' + Title.Caption = 'Value' Title.LabelBrush.Style = bsClear Title.TextFormat = tfHTML end> @@ -218,7 +218,7 @@ object DistribFrm: TDistribFrm Top = 8 Width = 144 Caption = 'Normal Distribution' - OnClick = NDChkClick + OnClick = DistributionClick TabOrder = 0 end object tChk: TRadioButton @@ -227,9 +227,8 @@ object DistribFrm: TDistribFrm Top = 35 Width = 144 Caption = 'Student t Distribution' - OnClick = tChkClick + OnClick = DistributionClick TabOrder = 3 - Visible = False end object FChk: TRadioButton Left = 16 @@ -237,7 +236,7 @@ object DistribFrm: TDistribFrm Top = 62 Width = 144 Caption = 'Central F Distribution' - OnClick = FChkClick + OnClick = DistributionClick TabOrder = 2 end object ChiChk: TRadioButton @@ -246,7 +245,7 @@ object DistribFrm: TDistribFrm Top = 89 Width = 144 Caption = 'Chi-Square Distribution' - OnClick = ChiChkClick + OnClick = DistributionClick TabOrder = 1 end end diff --git a/applications/lazstats/source/forms/simulations/distribunit.pas b/applications/lazstats/source/forms/simulations/distribunit.pas index 2cdd96ea2..c70f56be6 100644 --- a/applications/lazstats/source/forms/simulations/distribunit.pas +++ b/applications/lazstats/source/forms/simulations/distribunit.pas @@ -51,25 +51,22 @@ type DF2Label: TLabel; MeanLabel: TLabel; GroupBox1: TGroupBox; - procedure ChiChkClick(Sender: TObject); procedure ComputeBtnClick(Sender: TObject); - procedure FChkClick(Sender: TObject); procedure FormActivate(Sender: TObject); procedure FormShow(Sender: TObject); procedure CalcChiSq(const AX: Double; out AY: Double); procedure CalcF(const AX: Double; out AY: Double); procedure CalcND(const AX: Double; out AY: Double); procedure Calct(const AX: Double; out AY: Double); - procedure NDChkClick(Sender: TObject); + procedure DistributionClick(Sender: TObject); procedure PrintBtnClick(Sender: TObject); procedure ResetBtnClick(Sender: TObject); procedure SaveBtnClick(Sender: TObject); - procedure tChkClick(Sender: TObject); private { private declarations } DF1: Integer; DF2: Integer; - procedure NDPlot; + procedure NormalDistPlot; procedure ChiPlot; procedure FPlot; procedure tPlot; @@ -85,9 +82,8 @@ var implementation uses - //spe, // a numlib unit (tDist) - OSPrinters, - TAChartUtils, TADrawerSVG, TAPrint; + TAChartUtils, TADrawerSVG, TAPrint, + MathUnit; { TDistribFrm } @@ -131,28 +127,8 @@ begin end; procedure TDistribFrm.CalcF(const AX: Double; out AY: Double); -var - ratio1, ratio2, ratio3, ratio4: double; - part1, part2, part3, part4, part5, part6, part7, part8, part9: double; begin - // Returns the height of the density curve for the F statistic - ratio1 := (DF1 + DF2) / 2.0; - ratio2 := (DF1 - 2.0) / 2.0; - ratio3 := DF1 / 2.0; - ratio4 := DF2 / 2.0; - part1 := exp(lngamma(ratio1)); - part2 := power(DF1, ratio3); - part3 := power(DF2, ratio4); - part4 := exp(lngamma(ratio3)); - part5 := exp(lngamma(ratio4)); - part6 := power(AX, ratio2); - part7 := power((AX*DF1 + DF2), ratio1); - part8 := (part1 * part2 * part3) / (part4 * part5); - if (part7 = 0.0) then - part9 := 0.0 - else - part9 := part6 / part7; - AY := part8 * part9; + AY := FDensity(AX, DF1, DF2); end; procedure TDistribFrm.CalcND(const AX: Double; out AY: Double); @@ -162,17 +138,35 @@ end; procedure TDistribFrm.CalcChiSq(const AX: Double; out AY: Double); begin - AY := 1.0 / (power(2.0, DF1*0.5) * exp(lngamma(DF1*0.5))) * power(AX, (DF1-2.0)*0.5) * (1.0 / exp(AX*0.5)); + AY := Chi2Density(AX, DF1); end; procedure TDistribFrm.Calct(const AX: Double; out AY: Double); begin - AY := Student(AX, DF1, 1.0); - //AY := tDist(AX, DF1, 1); + AY := tDensity(AX, DF1); end; -procedure TDistribFrm.NDChkClick(Sender: TObject); +procedure TDistribFrm.DistributionClick(Sender: TObject); +var + rb: TRadiobutton; begin + rb := Sender as TRadioButton; + + GroupBox2.Enabled := rb.Checked; + + AlphaLabel.Enabled := rb.Checked; + AlphaEdit.Enabled := rb.Checked; + + DF1Edit.Enabled := (rb <> NDChk) and rb.Checked; + DF1Label.Enabled := DF1Edit.Enabled; + + DF2Edit.Enabled := (rb = FChk) and rb.Checked; + DF2Label.Enabled := DF2Edit.Enabled; + + MeanEdit.Enabled := false; + MeanLabel.Enabled := false; + { + if NDChk.Checked then begin GroupBox2.Enabled := true; @@ -187,6 +181,7 @@ begin end else GroupBox2.Enabled := false; + } end; procedure TDistribFrm.PrintBtnClick(Sender: TObject); @@ -233,7 +228,7 @@ begin ok := false; if NDChk.Checked then begin - NDPlot(); + NormalDistPlot(); ok := true; end; @@ -259,60 +254,7 @@ begin MessageDlg('Please select a distribution.', mtError, [mbOK], 0); end; -procedure TDistribFrm.FChkClick(Sender: TObject); -begin - if FChk.Checked then - begin - GroupBox2.Enabled := true; - DF2Label.Enabled := true; - AlphaLabel.Enabled := true; - AlphaEdit.Enabled := true; - DF1Edit.Enabled := true; - DF2Edit.Enabled := true; - DF1Label.Enabled := true; - MeanLabel.Enabled := false; - MeanEdit.Enabled := false; - end - else - GroupBox2.Enabled := false; -end; - -procedure TDistribFrm.ChiChkClick(Sender: TObject); -begin - if ChiChk.Checked then - begin - GroupBox2.Enabled := true; - DF1Label.Enabled := true; - DF1Edit.Enabled := true; - DF2Label.Enabled := false; - MeanLabel.Enabled := false; - AlphaLabel.Enabled := true; - AlphaEdit.Enabled := true; - DF2Edit.Enabled := false; - MeanEdit.Enabled := false; - end else - GroupBox2.Enabled := false; -end; - -procedure TDistribFrm.tChkClick(Sender: TObject); -begin - if tChk.Checked then - begin - GroupBox2.Enabled := true; - DF1Label.Enabled := true; - DF1Edit.Enabled := true; - DF2Label.Enabled := false; - MeanLabel.Enabled := false; - AlphaLabel.Enabled := true; - AlphaEdit.Enabled := true; - DF2Edit.Enabled := false; - MeanEdit.Enabled := false; - end else - GroupBox2.Enabled := false; -end; - - -procedure TDistribFrm.NDPlot; +procedure TDistribFrm.NormalDistPlot; var alpha: Double; zMax, zCrit, pCrit: Double; @@ -416,7 +358,7 @@ begin Chart.Title.Text.Add(Format('Critical value = %.3f', [tCrit])); Chart.Title.Visible := true; Chart.BottomAxis.Title.Caption := 't'; - FuncSeries.Extent.XMin := -tMax; + FuncSeries.Extent.XMin := -tmax; FuncSeries.Extent.XMax := tMax; FuncSeries.Extent.UseXMin := true; FuncSeries.Extent.UseXMax := true; diff --git a/applications/lazstats/source/units/functionslib.pas b/applications/lazstats/source/units/functionslib.pas index 45fac02af..c2b28c56a 100644 --- a/applications/lazstats/source/units/functionslib.pas +++ b/applications/lazstats/source/units/functionslib.pas @@ -10,7 +10,6 @@ uses MainUnit, dataprocs; function chisquaredprob(X : double; k : integer) : double; -function gammln(xx : double) : double; PROCEDURE matinv(VAR a, vtimesw, v, w: DblDyneMat; n: integer); FUNCTION sign(a,b: double): double; FUNCTION isign(a,b : integer): integer; @@ -66,6 +65,9 @@ procedure poisson_pdf ( x : integer; VAR a : double; VAR pdf : double ); implementation +uses + MathUnit; + function chisquaredprob(X : double; k : integer) : double; var factor : double; // factor which multiplies sum of series @@ -86,7 +88,7 @@ begin begin x1 := 0.5 * X; k1 := 0.5 * k; - g := gammln(k1 + 1); + g := GammaLn(k1 + 1); factor := exp(k1 * ln(x1) - g - x1); sum := 0.0; if factor > 0 then @@ -104,34 +106,7 @@ begin end; //end if .. else Result := chi2prob; end; -//--------------------------------------------------------------------- -function gammln(xx : double) : double; -var - X, tmp, ser : double; - cof : array[0..5] of double; - j : integer; - -begin - cof[0] := 76.18009173; - cof[1] := -86.50532033; - cof[2] := 24.01409822; - cof[3] := -1.231739516; - cof[4] := 0.00120858003; - cof[5] := -0.00000536382; - - X := xx - 1.0; - tmp := X + 5.5; - tmp := tmp - ((X + 0.5) * ln(tmp)); - ser := 1.0; - for j := 0 to 5 do - begin - X := X + 1.0; - ser := ser + cof[j] / X; - end; - Result := ( -tmp + ln(2.50662827465 * ser) ); -end; -//------------------------------------------------------------------- PROCEDURE matinv(VAR a, vtimesw, v, w: DblDyneMat; n: integer); (* adapted from the singular value decomposition of a matrix *) @@ -568,74 +543,6 @@ BEGIN END; //----------------------------------------------------------------- -function betacf(a,b,x: double): extended; -CONST - itmax=100; - eps=3.0e-7; -VAR - tem,qap,qam,qab,em,d: extended; - bz,bpp,bp,bm,az,app: extended; - am,aold,ap: extended; - m: integer; -BEGIN - am := 1.0; - bm := 1.0; - az := 1.0; - qab := a+b; - qap := a+1.0; - qam := a-1.0; - bz := 1.0 - qab * x / qap; - FOR m := 1 to itmax DO BEGIN - em := m; - tem := em+em; - d := em * (b - m) * x / ((qam + tem) * (a + tem)); - ap := az + d * am; - bp := bz + d * bm; - term1 := -(a + em); - term2 := qab + em; - term3 := term1 * term2 * x; - term4 := a + tem; - term5 := qap + tem; - term6 := term4 * term5; - d := term3 / term6; - app := ap + d * az; - bpp := bp + d * bz; - aold := az; - am := ap/bpp; - bm := bp/bpp; - az := app/bpp; - bz := 1.0; - IF ((abs(az-aold)) < (eps*abs(az))) THEN - Break; - END; - { ShowMessage('WARNING! a or b too big, or itmax too small in betacf');} - Result := az -END; - -FUNCTION betai(a,b,x: extended): extended; -VAR - bt: extended; -BEGIN - IF ((x <= 0.0) OR (x >= 1.0)) THEN BEGIN - { ShowMessage('ERROR! Problem in routine BETAI');} - betai := 0.5; - exit; - END; - IF ((x <= 0.0) OR (x >= 1.0)) THEN bt := 0.0 - ELSE - begin - term1 := gammln(a + b) - - gammln(a) - gammln(b); - term2 := a * ln(x); - term3 := b * ln(1.0 - x); - term4 := term1 + term2 + term3; - bt := exp(term4); - term5 := (a + 1.0) / (a + b + 2.0); - end; - IF x < term5 then betai := bt * betacf(a,b,x) / a - ELSE betai := 1.0 - bt * betacf(b,a,1.0-x) / b -END; - begin { fprob function } if f <= 0.0 then probf := 1.0 else probf := (betai(0.5*df2,0.5*df1,df2/(df2+df1*f)) + diff --git a/applications/lazstats/source/units/mathunit.pas b/applications/lazstats/source/units/mathunit.pas new file mode 100644 index 000000000..5a16c0812 --- /dev/null +++ b/applications/lazstats/source/units/mathunit.pas @@ -0,0 +1,188 @@ +unit MathUnit; + +{ extract some math functions from functionslib for easier testing } + +{$mode objfpc}{$H+} + +interface + +uses + Classes, SysUtils; + +function Beta(a, b: Double): Extended; +function BetaI(a,b,x: Double): Extended; + +function GammaLn(x: double): Extended; + +function tDist(x: Double; N: Integer; OneSided: Boolean): Double; +function tDensity(x: Double; N: Integer): Double; + +function FDensity(x: Double; DF1, DF2: Integer): Double; + +function Chi2Density(x: Double; N: Integer): Double; + +implementation + +uses + Math; + +function Beta(a, b: Double): Extended; +begin + if (a > 0) and (b > 0) then + Result := exp(GammaLn(a) + GammaLn(b) - GammaLn(a+b)) + else + raise Exception.Create('Invalid argument for beta function.'); +end; + +function BetaCF(a, b, x: double): Extended; +const + itmax = 100; + eps = 3.0e-7; +var + tem,qap,qam,qab,em,d: extended; + bz,bpp,bp,bm,az,app: extended; + am,aold,ap: extended; + term1, term2, term3, term4, term5, term6: extended; + m: integer; +BEGIN + am := 1.0; + bm := 1.0; + az := 1.0; + qab := a+b; + qap := a+1.0; + qam := a-1.0; + bz := 1.0 - qab * x / qap; + for m := 1 to itmax do begin + em := m; + tem := em+em; + d := em * (b - m) * x / ((qam + tem) * (a + tem)); + ap := az + d * am; + bp := bz + d * bm; + term1 := -(a + em); + term2 := qab + em; + term3 := term1 * term2 * x; + term4 := a + tem; + term5 := qap + tem; + term6 := term4 * term5; + d := term3 / term6; + app := ap + d * az; + bpp := bp + d * bz; + aold := az; + am := ap/bpp; + bm := bp/bpp; + az := app/bpp; + bz := 1.0; + if ((abs(az-aold)) < (eps*abs(az))) then + Break; + end; + + { ShowMessage('WARNING! a or b too big, or itmax too small in betacf');} + Result := az +end; + +function BetaI(a,b,x: Double): extended; +var + bt: extended; + term1, term2, term3, term4, term5: extended; +begin + if ((x < 0.0) or (x > 1.0)) then begin + { ShowMessage('ERROR! Problem in routine BETAI');} + Result := 0.5; + exit; + end; + + if ((x <= 0.0) or (x >= 1.0)) then + bt := 0.0 + else + begin + term1 := GammaLn(a + b) - GammaLn(a) - GammaLn(b); + term2 := a * ln(x); + term3 := b * ln(1.0 - x); + term4 := term1 + term2 + term3; + bt := exp(term4); + end; + term5 := (a + 1.0) / (a + b + 2.0); + if x < term5 then + Result := bt * betacf(a,b,x) / a + else + Result := 1.0 - bt * betacf(b,a,1.0-x) / b +end; + +function GammaLn(x: double): Extended; +var + tmp, ser: double; + cof: array[0..5] of double; + j: integer; +begin + cof[0] := 76.18009173; + cof[1] := -86.50532033; + cof[2] := 24.01409822; + cof[3] := -1.231739516; + cof[4] := 0.00120858003; + cof[5] := -0.00000536382; + + x := x - 1.0; + tmp := x + 5.5; + tmp := tmp - (x + 0.5) * ln(tmp); + ser := 1.0; + for j := 0 to 5 do + begin + x := x + 1.0; + ser := ser + cof[j] / x; + end; + Result := -tmp + ln(2.50662827465 * ser); +end; + +// Calculates the (cumulative) t distribution function for N degrees of freedom +function tDist(x: Double; N: Integer; OneSided: Boolean): Double; +begin + Result := betai(0.5*N, 0.5, N/(N + sqr(x))); + if OneSided then Result := Result * 0.5; +end; + +// Returns the density curve for the t statistic with N degrees of freedom +function tDensity(x: Double; N: Integer): Double; +var + factor: Double; +begin + factor := exp(gammaLn((N+1)/2) - gammaLn(N/2)) / sqrt(N * pi); + Result := factor * Power(1 + sqr(x)/N, (1-n)/2); +end; + +// Returns the density curve for the F statistic for DF1 and DF2 degrees of freedom. +function FDensity(x: Double; DF1, DF2: Integer): Double; +var + ratio1, ratio2, ratio3, ratio4: double; + part1, part2, part3, part4, part5, part6, part7, part8, part9: double; +begin + ratio1 := (DF1 + DF2) / 2.0; + ratio2 := (DF1 - 2.0) / 2.0; + ratio3 := DF1 / 2.0; + ratio4 := DF2 / 2.0; + part1 := exp(gammaln(ratio1)); + part2 := power(DF1, ratio3); + part3 := power(DF2, ratio4); + part4 := exp(gammaln(ratio3)); + part5 := exp(gammaln(ratio4)); + part6 := power(x, ratio2); + part7 := power((x*DF1 + DF2), ratio1); + part8 := (part1 * part2 * part3) / (part4 * part5); + if (part7 = 0.0) then + part9 := 0.0 + else + part9 := part6 / part7; + Result := part8 * part9; +end; + + +// Returns the density curve of the chi2 statistic for N degrees of freedom +function Chi2Density(x: Double; N: Integer): Double; +var + factor: Double; +begin + factor := Power(2.0, N * 0.5) * exp(gammaLN(N * 0.5)); + Result := power(x, (N-2.0) * 0.5) / (factor * exp(x * 0.5)); +end; + +end. +