diff --git a/applications/lazstats/source/LazStats.lpi b/applications/lazstats/source/LazStats.lpi index 07a93c1d3..56254e1e6 100644 --- a/applications/lazstats/source/LazStats.lpi +++ b/applications/lazstats/source/LazStats.lpi @@ -1429,6 +1429,11 @@ + + + + + diff --git a/applications/lazstats/source/forms/simulations/distribunit.lfm b/applications/lazstats/source/forms/simulations/distribunit.lfm index 7e49c860e..312df149a 100644 --- a/applications/lazstats/source/forms/simulations/distribunit.lfm +++ b/applications/lazstats/source/forms/simulations/distribunit.lfm @@ -2,12 +2,12 @@ object DistribFrm: TDistribFrm Left = 420 Height = 432 Top = 215 - Width = 687 + Width = 701 HelpType = htKeyword HelpKeyword = 'html/DistributionPlotsandCriticalValu.htm' Caption = 'Distributions' ClientHeight = 432 - ClientWidth = 687 + ClientWidth = 701 OnActivate = FormActivate OnShow = FormShow Position = poMainFormCenter @@ -19,7 +19,7 @@ object DistribFrm: TDistribFrm AnchorSideRight.Side = asrBottom AnchorSideBottom.Control = Owner AnchorSideBottom.Side = asrBottom - Left = 624 + Left = 638 Height = 25 Top = 401 Width = 55 @@ -35,7 +35,7 @@ object DistribFrm: TDistribFrm object ComputeBtn: TButton AnchorSideTop.Control = CloseBtn AnchorSideRight.Control = CloseBtn - Left = 540 + Left = 554 Height = 25 Top = 401 Width = 76 @@ -50,7 +50,7 @@ object DistribFrm: TDistribFrm object ResetBtn: TButton AnchorSideTop.Control = CloseBtn AnchorSideRight.Control = ComputeBtn - Left = 478 + Left = 492 Height = 25 Top = 401 Width = 54 @@ -72,7 +72,7 @@ object DistribFrm: TDistribFrm Left = 0 Height = 8 Top = 385 - Width = 687 + Width = 701 Anchors = [akLeft, akRight, akBottom] Shape = bsBottomLine end @@ -84,7 +84,7 @@ object DistribFrm: TDistribFrm Left = 8 Height = 377 Top = 8 - Width = 479 + Width = 454 Anchors = [akTop, akLeft, akRight, akBottom] BorderSpacing.Left = 8 BorderSpacing.Top = 8 @@ -92,7 +92,7 @@ object DistribFrm: TDistribFrm BevelOuter = bvNone BorderStyle = bsSingle ClientHeight = 373 - ClientWidth = 475 + ClientWidth = 450 Color = clWhite ParentColor = False TabOrder = 3 @@ -100,7 +100,7 @@ object DistribFrm: TDistribFrm Left = 6 Height = 361 Top = 6 - Width = 463 + Width = 438 AxisList = < item Grid.Color = clSilver @@ -146,6 +146,12 @@ object DistribFrm: TDistribFrm end object VertLineSeries: TLineSeries Active = False + LinePen.Color = clRed + LinePen.Style = psDot + end + object HorLineSeries: TLineSeries + LinePen.Color = clRed + LinePen.Style = psDot end end end @@ -181,42 +187,42 @@ object DistribFrm: TDistribFrm AnchorSideRight.Control = Owner AnchorSideRight.Side = asrBottom AnchorSideBottom.Control = Bevel1 - Left = 499 + Left = 474 Height = 377 Top = 8 - Width = 180 + Width = 219 Anchors = [akTop, akRight, akBottom] AutoSize = True BorderSpacing.Top = 8 BorderSpacing.Right = 8 BevelOuter = bvNone ClientHeight = 377 - ClientWidth = 180 + ClientWidth = 219 TabOrder = 6 object GroupBox1: TGroupBox AnchorSideLeft.Control = ParameterPanel AnchorSideTop.Control = ParameterPanel AnchorSideRight.Side = asrBottom Left = 0 - Height = 136 + Height = 175 Top = 0 - Width = 180 + Width = 219 AutoSize = True - Caption = 'Plot Distribution:' + Caption = 'Plot Distribution' ChildSizing.LeftRightSpacing = 16 ChildSizing.TopBottomSpacing = 8 ChildSizing.VerticalSpacing = 8 - ChildSizing.EnlargeVertical = crsHomogenousSpaceResize + ChildSizing.EnlargeVertical = crsScaleChilds ChildSizing.Layout = cclLeftToRightThenTopToBottom ChildSizing.ControlsPerLine = 1 - ClientHeight = 116 - ClientWidth = 176 + ClientHeight = 155 + ClientWidth = 215 TabOrder = 0 object NDChk: TRadioButton Left = 16 Height = 19 Top = 8 - Width = 144 + Width = 183 Caption = 'Normal Distribution' OnClick = DistributionClick TabOrder = 0 @@ -225,8 +231,8 @@ object DistribFrm: TDistribFrm Left = 16 Height = 19 Top = 35 - Width = 144 - Caption = 'Student t Distribution' + Width = 183 + Caption = 'Student t Distribution (1-sided)' OnClick = DistributionClick TabOrder = 3 end @@ -234,7 +240,7 @@ object DistribFrm: TDistribFrm Left = 16 Height = 19 Top = 62 - Width = 144 + Width = 183 Caption = 'Central F Distribution' OnClick = DistributionClick TabOrder = 2 @@ -243,11 +249,28 @@ object DistribFrm: TDistribFrm Left = 16 Height = 19 Top = 89 - Width = 144 + Width = 183 Caption = 'Chi-Square Distribution' OnClick = DistributionClick TabOrder = 1 end + object Bevel2: TBevel + Left = 16 + Height = 4 + Top = 116 + Width = 183 + Constraints.MaxHeight = 4 + Shape = bsBottomLine + end + object CumulativeChk: TCheckBox + Left = 16 + Height = 19 + Top = 128 + Width = 183 + Caption = 'Cumulative' + OnChange = CumulativeChkChange + TabOrder = 4 + end end object GroupBox2: TGroupBox AnchorSideLeft.Control = GroupBox1 @@ -256,21 +279,22 @@ object DistribFrm: TDistribFrm AnchorSideRight.Control = GroupBox1 AnchorSideRight.Side = asrBottom Left = 0 - Height = 140 - Top = 152 - Width = 180 + Height = 107 + Top = 191 + Width = 219 Anchors = [akTop, akLeft, akRight] + AutoSize = True BorderSpacing.Top = 16 BorderSpacing.Bottom = 12 Caption = 'Parameters' - ClientHeight = 120 - ClientWidth = 176 + ClientHeight = 87 + ClientWidth = 215 TabOrder = 1 object AlphaLabel: TLabel AnchorSideTop.Control = AlphaEdit AnchorSideTop.Side = asrCenter AnchorSideRight.Control = AlphaEdit - Left = 29 + Left = 68 Height = 15 Top = 6 Width = 84 @@ -283,7 +307,7 @@ object DistribFrm: TDistribFrm AnchorSideTop.Control = DF1Edit AnchorSideTop.Side = asrCenter AnchorSideRight.Control = DF1Edit - Left = 22 + Left = 61 Height = 15 Top = 33 Width = 91 @@ -297,7 +321,7 @@ object DistribFrm: TDistribFrm AnchorSideTop.Control = DF2Edit AnchorSideTop.Side = asrCenter AnchorSideRight.Control = DF2Edit - Left = 22 + Left = 61 Height = 15 Top = 60 Width = 91 @@ -307,25 +331,11 @@ object DistribFrm: TDistribFrm Caption = 'Deg. Freedom (2)' ParentColor = False end - object MeanLabel: TLabel - AnchorSideTop.Control = MeanEdit - AnchorSideTop.Side = asrCenter - AnchorSideRight.Control = MeanEdit - Left = 83 - Height = 15 - Top = 93 - Width = 30 - Anchors = [akTop, akRight] - BorderSpacing.Right = 8 - Caption = 'Mean' - ParentColor = False - Visible = False - end object AlphaEdit: TEdit AnchorSideTop.Control = GroupBox2 AnchorSideRight.Control = GroupBox2 AnchorSideRight.Side = asrBottom - Left = 121 + Left = 160 Height = 23 Top = 2 Width = 43 @@ -341,7 +351,7 @@ object DistribFrm: TDistribFrm AnchorSideTop.Side = asrBottom AnchorSideRight.Control = AlphaEdit AnchorSideRight.Side = asrBottom - Left = 121 + Left = 160 Height = 23 Top = 29 Width = 43 @@ -351,40 +361,36 @@ object DistribFrm: TDistribFrm TabOrder = 1 Text = 'DF1Edit' end - object MeanEdit: TEdit - AnchorSideTop.Control = DF2Edit - AnchorSideTop.Side = asrBottom - AnchorSideRight.Control = AlphaEdit - AnchorSideRight.Side = asrBottom - Left = 121 - Height = 23 - Top = 89 - Width = 43 - Alignment = taRightJustify - Anchors = [akTop, akRight] - BorderSpacing.Top = 4 - BorderSpacing.Bottom = 8 - TabOrder = 3 - Text = 'MeanEdit' - Visible = False - end object DF2Edit: TEdit AnchorSideTop.Control = DF1Edit AnchorSideTop.Side = asrBottom AnchorSideRight.Control = AlphaEdit AnchorSideRight.Side = asrBottom - Left = 121 + Left = 160 Height = 23 Top = 56 Width = 43 Alignment = taRightJustify Anchors = [akTop, akRight] BorderSpacing.Top = 4 - BorderSpacing.Bottom = 10 + BorderSpacing.Bottom = 8 TabOrder = 2 Text = 'DF2Edit' end end + object ShowCriticalValuesChk: TCheckBox + AnchorSideLeft.Control = ParameterPanel + AnchorSideTop.Control = GroupBox2 + AnchorSideTop.Side = asrBottom + Left = 0 + Height = 19 + Top = 310 + Width = 123 + Caption = 'Show critical values' + Checked = True + State = cbChecked + TabOrder = 2 + end end object SavePictureDialog: TSavePictureDialog Filter = 'Graphic (*.png;*.bmp;*.jpeg;*.jpg;*.jpe;*.jfif;*.svg)|*.png;*.bmp;*.jpeg;*.jpg;*.jpe;*.jfif;*.svg|Portable Network Graphic (*.png)|*.png|Bitmaps (*.bmp)|*.bmp|Joint Picture Expert Group (*.jpeg;*.jpg;*.jpe;*.jfif)|*.jpeg;*.jpg;*.jpe;*.jfif|Scaleable Vector Graphic (*.svg)|*.svg|All Files (*.*)|*.*' diff --git a/applications/lazstats/source/forms/simulations/distribunit.pas b/applications/lazstats/source/forms/simulations/distribunit.pas index c70f56be6..fb9b341ad 100644 --- a/applications/lazstats/source/forms/simulations/distribunit.pas +++ b/applications/lazstats/source/forms/simulations/distribunit.pas @@ -23,7 +23,11 @@ type TDistribFrm = class(TForm) AlphaEdit: TEdit; Bevel1: TBevel; + Bevel2: TBevel; + ShowCriticalValuesChk: TCheckBox; + HorLineSeries: TLineSeries; ChartToolset: TChartToolset; + CumulativeChk: TCheckBox; PanDragTool: TPanDragTool; tChk: TRadioButton; ZoomDragTool: TZoomDragTool; @@ -39,7 +43,6 @@ type DF1Edit: TEdit; DF2Edit: TEdit; FChk: TRadioButton; - MeanEdit: TEdit; NDChk: TRadioButton; ChartPanel: TPanel; ResetBtn: TButton; @@ -49,15 +52,19 @@ type AlphaLabel: TLabel; DF1Label: TLabel; DF2Label: TLabel; - MeanLabel: TLabel; GroupBox1: TGroupBox; procedure ComputeBtnClick(Sender: TObject); + procedure CumulativeChkChange(Sender: TObject); procedure FormActivate(Sender: TObject); procedure FormShow(Sender: TObject); - procedure CalcChiSq(const AX: Double; out AY: Double); + procedure CalcChi2(const AX: Double; out AY: Double); + procedure CalcChi2_Cumulative(const AX: Double; out AY: Double); procedure CalcF(const AX: Double; out AY: Double); + procedure CalcF_Cumulative(const AX: Double; out AY: Double); procedure CalcND(const AX: Double; out AY: Double); - procedure Calct(const AX: Double; out AY: Double); + procedure CalcND_Cumulative(const AX: Double; out AY: Double); + procedure CalcT(const AX: Double; out AY: Double); + procedure CalcT_Cumulative(const AX: Double; out AY: Double); procedure DistributionClick(Sender: TObject); procedure PrintBtnClick(Sender: TObject); procedure ResetBtnClick(Sender: TObject); @@ -67,7 +74,7 @@ type DF1: Integer; DF2: Integer; procedure NormalDistPlot; - procedure ChiPlot; + procedure Chi2Plot; procedure FPlot; procedure tPlot; @@ -97,7 +104,6 @@ begin AlphaEdit.Text := FormatFloat('0.00', DEFAULT_ALPHA_LEVEL); DF1Edit.Text := ''; DF2Edit.Text := ''; - MeanEdit.Text := ''; GroupBox2.Enabled := false; FuncSeries.OnCalculate := nil; VertLineSeries.Active := false; @@ -105,6 +111,7 @@ begin Chart.BottomAxis.Title.Caption := 'Scale'; end; + procedure TDistribFrm.SaveBtnClick(Sender: TObject); var ext: String; @@ -121,31 +128,67 @@ begin end; end; + procedure TDistribFrm.FormShow(Sender: TObject); begin ResetBtnClick(self); end; + procedure TDistribFrm.CalcF(const AX: Double; out AY: Double); begin AY := FDensity(AX, DF1, DF2); end; + +procedure TDistribFrm.CalcF_Cumulative(const AX: Double; out AY: Double); +begin + AY := 1.0 - ProbF(AX, DF1, DF2); +end; + + procedure TDistribFrm.CalcND(const AX: Double; out AY: Double); begin AY := 1.0 / sqrt(TWO_PI) * exp(-sqr(AX)/ 2.0); end; -procedure TDistribFrm.CalcChiSq(const AX: Double; out AY: Double); + +procedure TDistribFrm.CalcND_Cumulative(const AX: Double; out AY: Double); +begin + // AY := ProbZ(AX); -- very slow + AY := NormalDist(AX); // borrowed from NumLib +end; + + +procedure TDistribFrm.CalcChi2(const AX: Double; out AY: Double); begin AY := Chi2Density(AX, DF1); end; + +procedure TDistribFrm.CalcChi2_Cumulative(const AX: Double; out AY: Double); +begin + AY := ChiSquaredProb(AX, DF1); +end; + + procedure TDistribFrm.Calct(const AX: Double; out AY: Double); begin AY := tDensity(AX, DF1); end; + +procedure TDistribFrm.CalcT_Cumulative(const AX: Double; out AY: Double); +const + ONE_SIDED = true; +begin + if AX < 0 then + AY := tDist(-AX, DF1, ONE_SIDED) + else + AY := 1.0 - tDist(AX, DF1, ONE_SIDED); +end; + + procedure TDistribFrm.DistributionClick(Sender: TObject); var rb: TRadiobutton; @@ -162,28 +205,9 @@ begin 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; - AlphaLabel.Enabled := true; - AlphaEdit.Enabled := true; - DF1Edit.Enabled := false; - DF1Label.Enabled := false; - DF2Label.Enabled := false; - MeanLabel.Enabled := false; - DF2Edit.Enabled := false; - MeanEdit.Enabled := false; - end - else - GroupBox2.Enabled := false; - } end; + procedure TDistribFrm.PrintBtnClick(Sender: TObject); const MARGIN = 10; @@ -212,6 +236,7 @@ begin end; end; + procedure TDistribFrm.ComputeBtnClick(Sender: TObject); var msg: String; @@ -240,7 +265,7 @@ begin if ChiChk.Checked then begin - ChiPlot(); + Chi2Plot(); ok := true; end; @@ -254,13 +279,23 @@ begin MessageDlg('Please select a distribution.', mtError, [mbOK], 0); end; +procedure TDistribFrm.CumulativeChkChange(Sender: TObject); +begin + if CumulativeChk.Checked then + Chart.LeftAxis.Title.Caption := 'Cumulative probability' + else + Chart.LeftAxis.Title.Caption := 'Probability density'; +end; + + procedure TDistribFrm.NormalDistPlot; var alpha: Double; - zMax, zCrit, pCrit: Double; + zMax, zMin, zCrit, pCrit: Double; begin alpha := StrToFloat(AlphaEdit.Text); zMax := InverseZ(0.9999); + zMin := -zMax; zCrit := inversez(1.0 - alpha); CalcND(zCrit, pCrit); @@ -270,28 +305,44 @@ begin Chart.Title.Text.Add(Format('Critical value = %.3f', [zCrit])); Chart.Title.Visible := true; Chart.BottomAxis.Title.Caption := 'z'; - FuncSeries.Extent.XMin := -zMax; - FuncSeries.Extent.XMax := +zMax; + FuncSeries.Extent.XMin := zMin; + FuncSeries.Extent.XMax := zMax; FuncSeries.Extent.UseXMin := true; FuncSeries.Extent.UseXMax := true; - FuncSeries.OnCalculate := @CalcND; + if CumulativeChk.Checked then + FuncSeries.OnCalculate := @CalcND_Cumulative + else + FuncSeries.OnCalculate := @CalcND; FuncSeries.DomainExclusions.Clear; - VertLineSeries.Clear; - VertLineSeries.AddXY(zCrit, 0); - VertLineSeries.AddXY(zCrit, pCrit); - VertLineSeries.Active := true; + + VertlineSeries.Clear; + if CumulativeChk.Checked then + begin + HorLineSeries.Clear; + HorLineSeries.AddXY(zMin, 1.0 - alpha); + HorLineSeries.AddXY(zCrit, 1.0 - alpha); + VertlineSeries.AddXY(zCrit, 1.0 - alpha); + VertLineSeries.AddXY(zCrit, 0); + end else + begin + VertLineSeries.AddXY(zCrit, 0); + VertLineSeries.AddXY(zCrit, pCrit); + end; + HorLineSeries.Active := ShowCriticalValuesChk.Checked and CumulativeChk.Checked; + VertLineSeries.Active := ShowCriticalValuesChk.Checked; end; -procedure TDistribFrm.ChiPlot; + +procedure TDistribFrm.Chi2Plot; var alpha: Double; - Chi2Max, Chi2Crit, pCrit: Double; + chi2Max, chi2Crit, pCrit: Double; begin alpha := StrToFloat(AlphaEdit.Text); DF1 := StrToInt(DF1Edit.Text); - Chi2Max := InverseChi(0.9999, DF1); - Chi2Crit := InverseChi(1.0 - alpha, DF1); - CalcChiSq(Chi2Crit, pCrit); + chi2Max := InverseChi(0.9999, DF1); + chi2Crit := InverseChi(1.0 - alpha, DF1); + CalcChi2(chi2Crit, pCrit); Chart.Title.Text.Clear; Chart.Title.Text.Add('Chi-Squared Distribution'); @@ -300,17 +351,32 @@ begin Chart.Title.Visible := true; Chart.BottomAxis.Title.Caption := 'χ2'; FuncSeries.Extent.XMin := 0; - FuncSeries.Extent.XMax := Chi2Max; + FuncSeries.Extent.XMax := chi2Max; FuncSeries.Extent.UseXMin := true; FuncSeries.Extent.UseXMax := true; - FuncSeries.OnCalculate := @CalcChiSq; + if CumulativeChk.Checked then + FuncSeries.OnCalculate := @CalcChi2_Cumulative + else + FuncSeries.OnCalculate := @CalcChi2; FuncSeries.DomainExclusions.AddRange(-Infinity, 0, [ioOpenEnd]); + VertLineSeries.Clear; - VertLineSeries.AddXY(Chi2Crit, 0); - VertLineSeries.AddXY(Chi2Crit, pCrit); - VertLineSeries.Active := true; + if CumulativeChk.Checked then begin + HorLineSeries.Clear; + HorLineSeries.AddXY(0, 1.0 - alpha); + HorLineSeries.AddXY(chi2Crit, 1.0 - alpha); + VertlineSeries.AddXY(chi2Crit, 1.0 - alpha); + VertLineSeries.AddXY(chi2Crit, 0); + end else + begin + VertLineSeries.AddXY(chi2Crit, 0); + VertLineSeries.AddXY(chi2Crit, pCrit); + end; + HorLineSeries.Active := ShowCriticalValuesChk.Checked and CumulativeChk.Checked; + VertLineSeries.Active := ShowCriticalValuesChk.Checked; end; + procedure TDistribFrm.FPlot; var alpha: Double; @@ -333,22 +399,39 @@ begin FuncSeries.Extent.XMax := FMax; FuncSeries.Extent.UseXMin := true; FuncSeries.Extent.UseXMax := true; - FuncSeries.OnCalculate := @CalcF; + if CumulativeChk.Checked then + FuncSeries.OnCalculate := @CalcF_Cumulative + else + FuncSeries.OnCalculate := @CalcF; FuncSeries.DomainExclusions.AddRange(-Infinity, 0, [ioOpenEnd]); + VertLineSeries.Clear; - VertLineSeries.AddXY(FCrit, 0); - VertLineSeries.AddXY(FCrit, pCrit); - VertLineSeries.Active := true; + if CumulativeChk.Checked then + begin + HorLineSeries.Clear; + HorLineSeries.AddXY(0, 1.0 - alpha); + HorLineSeries.AddXY(FCrit, 1.0 - alpha); + VertLineSeries.AddXY(FCrit, 1.0 - alpha); + VertLineSeries.AddXY(FCrit, 0); + end else + begin + VertLineSeries.AddXY(FCrit, 0); + VertLineSeries.AddXY(FCrit, pCrit); + end; + HorLineSeries.Active := ShowCriticalValuesChk.Checked and CumulativeChk.Checked; + VertLineSeries.Active := ShowCriticalValuesChk.Checked; end; + procedure TDistribFrm.tPlot; var alpha: Double; - tMax, tCrit, pCrit: Double; + tMin, tMax, tCrit, pCrit: Double; begin alpha := StrToFloat(AlphaEdit.Text); DF1 := StrToInt(DF1Edit.Text); tMax := Inverset(0.9999, DF1); + tMin := -tMax; tCrit := Inverset(1.0 - alpha, DF1); Calct(tCrit, pCrit); @@ -358,18 +441,34 @@ 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 := tMin; FuncSeries.Extent.XMax := tMax; FuncSeries.Extent.UseXMin := true; FuncSeries.Extent.UseXMax := true; - FuncSeries.OnCalculate := @Calct; + if CumulativeChk.Checked then + FuncSeries.OnCalculate := @CalcT_Cumulative + else + FuncSeries.OnCalculate := @CalcT; FuncSeries.DomainExclusions.Clear; + VertLineSeries.Clear; - VertLineSeries.AddXY(tCrit, 0); - VertLineSeries.AddXY(tCrit, pCrit); - VertLineSeries.Active := true; + if CumulativeChk.Checked then + begin + HorLineSeries.Clear; + HorLineSeries.AddXY(tMin, 1.0 - alpha); + HorLineSeries.AddXY(tCrit, 1.0 - alpha); + VertLineSeries.AddXY(tCrit, 1.0 - alpha); + VertLineSeries.AddXY(tCrit, 0); + end else + begin + VertLineSeries.AddXY(tCrit, 0); + VertLineSeries.AddXY(tCrit, pCrit); + end; + HorLineSeries.Active := ShowCriticalValuesChk.Checked and CumulativeChk.Checked; + VertLineSeries.Active := ShowCriticalValuesChk.Checked; end; + procedure TDistribFrm.FormActivate(Sender: TObject); var w: Integer; diff --git a/applications/lazstats/source/units/functionslib.pas b/applications/lazstats/source/units/functionslib.pas index c2b28c56a..079a4e1bb 100644 --- a/applications/lazstats/source/units/functionslib.pas +++ b/applications/lazstats/source/units/functionslib.pas @@ -980,16 +980,15 @@ begin end; //------------------------------------------------------------------- -function probt(t,df1 : double) : double; +// Returns the probability corresponding to a two-tailed t test. +function ProbT(t, df1: double): double; var - F, prob : double; + F, prob: double; begin - // Returns the probability corresponding to a two-tailed t test. - F := t * t; - prob := probf(F,1.0,df1); - Result := prob; + F := t * t; + prob := ProbF(F, 1.0, df1); + Result := prob; end; -//------------------------------------------------------------------------ function inverset(Probt, DF : double) : double; var diff --git a/applications/lazstats/source/units/mathunit.pas b/applications/lazstats/source/units/mathunit.pas index 5a16c0812..a0312a4d2 100644 --- a/applications/lazstats/source/units/mathunit.pas +++ b/applications/lazstats/source/units/mathunit.pas @@ -9,6 +9,10 @@ interface uses Classes, SysUtils; +function erf(x: Double): Double; +function erfc(x: Double) : Double; +function NormalDist(x: Double): Double; + function Beta(a, b: Double): Extended; function BetaI(a,b,x: Double): Extended; @@ -26,6 +30,104 @@ implementation uses Math; +// Calculates the error function +// /x +// erf(x) = 2/sqrt(pi) * | exp(-t²) dt +// /0 +// borrowed from NumLib +function erf(x: Double): Double; +const + xup = 6.25; + SQRT_PI = 1.7724538509055160; + c: array[1..18] of Double = ( + 1.9449071068178803e0, 4.20186582324414e-2, -1.86866103976769e-2, + 5.1281061839107e-3, -1.0683107461726e-3, 1.744737872522e-4, + -2.15642065714e-5, 1.7282657974e-6, -2.00479241e-8, + -1.64782105e-8, 2.0008475e-9, 2.57716e-11, + -3.06343e-11, 1.9158e-12, 3.703e-13, + -5.43e-14, -4.0e-15, 1.2e-15 + ); + d: array[1..17] of Double = ( + 1.4831105640848036e0, -3.010710733865950e-1, 6.89948306898316e-2, + -1.39162712647222e-2, 2.4207995224335e-3, -3.658639685849e-4, + 4.86209844323e-5, -5.7492565580e-6, 6.113243578e-7, + -5.89910153e-8, 5.2070091e-9, -4.232976e-10, + 3.18811e-11, -2.2361e-12, 1.467e-13, + -9.0e-15, 5.0e-16 + ); +var + t, s, s1, s2, x2: Double; + bovc, bovd, j: Integer; + sgn: Integer; +begin + bovc := SizeOf(c) div SizeOf(Double); + bovd := SizeOf(d) div SizeOf(Double); + t := abs(x); + if t <= 2 then + begin + x2 := sqr(x) - 2; + s1 := d[bovd]; + s2 := 0; + j := bovd - 1; + s := x2*s1 - s2 + d[j]; + while j > 1 do + begin + s2 := s1; + s1 := s; + j := j-1; + s := x2*s1 - s2 + d[j]; + end; + Result := (s - s2) * x / 2; + end else + if t < xup then + begin + x2 := 2 - 20 / (t+3); + s1 := c[bovc]; + s2 := 0; + j := bovc - 1; + s := x2*s1 - s2 + c[j]; + while j > 1 do + begin + s2 := s1; + s1 := s; + j := j-1; + s := x2*s1 - s2 + c[j]; + end; + x2 := ((s-s2) / (2*t)) * exp(-sqr(x)) / SQRT_PI; + if x < 0 then sgn := -1 else sgn := +1; + Result := (1 - x2) * sgn + end + else + if x < 0 then + Result := -1.0 + else + Result := +1.0; +end; + + +{ calculates the complementary error function erfc(x) = 1 - erf(x) } +function erfc(x: Double) : Double; +begin + Result := 1.0 - erf(x); +end; + + +// Cumulative normal distribution +// x = -INF ... INF --> 0 ... 1 +function NormalDist(x: Double): Double; +const + SQRT2 = sqrt(2.0); +begin + if x > 0 then + Result := (erf(x / SQRT2) + 1) * 0.5 + else + if x < 0 then + Result := (1.0 - erf(-x / SQRT2)) * 0.5 + else + Result := 0; +end; + + function Beta(a, b: Double): Extended; begin if (a > 0) and (b > 0) then