From e353889d67b797733992966859703f965d15a997 Mon Sep 17 00:00:00 2001 From: wp_xxyyzz Date: Sat, 26 Sep 2020 21:10:53 +0000 Subject: [PATCH] LazStats: Complete refactored NormalityUnit. git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7697 8e941d3f-bd1b-0410-a28a-d453659cc2b4 --- .../analysis/descriptive/normalityunit.lfm | 1 + .../analysis/descriptive/normalityunit.pas | 339 +++++++----------- 2 files changed, 127 insertions(+), 213 deletions(-) diff --git a/applications/lazstats/source/forms/analysis/descriptive/normalityunit.lfm b/applications/lazstats/source/forms/analysis/descriptive/normalityunit.lfm index 6e2c5e422..20e0ed072 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/normalityunit.lfm +++ b/applications/lazstats/source/forms/analysis/descriptive/normalityunit.lfm @@ -54,6 +54,7 @@ object NormalityFrm: TNormalityFrm BorderSpacing.Left = 8 Caption = 'Close' ModalResult = 11 + OnClick = CloseBtnClick TabOrder = 0 end object ComputeBtn: TButton diff --git a/applications/lazstats/source/forms/analysis/descriptive/normalityunit.pas b/applications/lazstats/source/forms/analysis/descriptive/normalityunit.pas index 6c3f0ab0c..bca6f8767 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/normalityunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/normalityunit.pas @@ -33,6 +33,7 @@ type VarOutBtn: TBitBtn; Label1: TLabel; VarList: TListBox; + procedure CloseBtnClick(Sender: TObject); procedure ComputeBtnClick(Sender: TObject); procedure FormActivate(Sender: TObject); procedure FormCreate(Sender: TObject); @@ -49,7 +50,6 @@ type function Calc_ShapiroWilks(const AData: DblDyneVec; out W, Prob: Double): Boolean; procedure Calc_Lilliefors(const AData: DblDyneVec; out ASkew, AKurtosis, AStat: Double; out AConclusion: String); - function Norm(z : double): double; procedure PlotData(AData: DblDyneVec); function PrepareData(const VarName: String): DblDyneVec; procedure UpdateBtnStates; @@ -67,84 +67,51 @@ implementation {$R *.lfm} uses - Math, + Math, MathUnit, TAChartUtils, TATextElements, TACustomSeries, TATransformations, - TACustomSource, TASources, + TAChartAxisUtils, TAChartAxis, TACustomSource, TASources, Utils; -{ TNormalityFrm } -procedure TNormalityFrm.FormActivate(Sender: TObject); -var - w: Integer; -begin - if FAutoSized then - exit; - - w := MaxValue([ResetBtn.Width, ComputeBtn.Width, CloseBtn.Width]); - ResetBtn.Constraints.MinWidth := w; - ComputeBtn.Constraints.MinWidth := w; - CloseBtn.Constraints.MinWidth := w; - - ParamsPanel.Constraints.MinWidth := Max( - 3*w + 2*CloseBtn.BorderSpacing.Left, - Max(Label1.Width, Label2.Width) + VarInBtn.Width + 2*VarList.BorderSpacing.Right); - ParamsPanel.Constraints.MinHeight := VarOutBtn.Top + VarOutBtn.Height + - VarOutBtn.BorderSpacing.Bottom + Bevel1.Height + Panel1.Height + - Panel1.BorderSpacing.Top; - - Constraints.MinWidth := ParamsPanel.Constraints.MinWidth + 300; - Constraints.MinHeight := ParamsPanel.Constraints.MinHeight + 2*ParamsPanel.BorderSpacing.Top; - - Position := poDesigned; - FAutoSized := True; -end; - -procedure TNormalityFrm.FormCreate(Sender: TObject); +function PopulateLeftMarks(AOwner: TComponent): TListChartSource; const - MARKS: array[0..23] of double = (0.001, 0.002, 0.005, 0.01, 0.02, 0.03, 0.05, 0.07, - 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.93, 0.95, 0.98, 0.99, - 0.995, 0.998, 0.999); + MARKS: array[0..12] of double = ( + 0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.4); var - T : TChartAxisTransformations; - InvNormDistTransform: TAxisTransform; - L: TListChartSource; i: Integer; begin - Assert(OS3MainFrm <> nil); - - InitForm(self); - - FReportFrame := TReportFrame.Create(self); - FReportFrame.Parent := ReportPage; - FReportFrame.Align := alClient; - - FChartFrame := TChartFrame.Create(self); - FChartFrame.Parent := ChartPage; - FChartFrame.Align := alClient; - FChartFrame.Chart.BottomAxis.Intervals.MaxLength := 80; - FChartFrame.Chart.BottomAxis.Intervals.MinLength := 30; - - T := TChartAxisTransformations.Create(FChartFrame); - InvNormDistTransform := TCumulNormDistrAxisTransform.Create(T); - InvNormDistTransform.Transformations := T; - FChartFrame.Chart.LeftAxis.Transformations := T; - - FChartFrame.Chart.LeftAxis.Intervals.Tolerance := 10; - FChartFrame.Chart.LeftAxis.Intervals.Count := 30; - FChartFrame.Chart.LeftAxis.Intervals.Options := FChartFrame.Chart.leftAxis.Intervals.options + [aipUseCount]; //aipGraphCoords]; - FChartFrame.Chart.LeftAxis.Marks.OverlapPolicy := opHideNeighbour; - { - L := TListChartSource.Create(FChartFrame.Chart); - for i := 0 to High(MARKS) do - L.Add(MARKS[i], MARKS[i]); - FChartFrame.Chart.LeftAxis.Marks.Source := L; - FChartFrame.Chart.LeftAxis.Marks.Style := smsLabel; - } - Reset; + Result := TListChartSource.Create(AOwner); + for i := 0 to High(MARKS) do Result.Add(MARKS[i], MARKS[i]); + Result.Add(0.5, 0.5); + for i := High(Marks) downto 0 do Result.Add(1-MARKS[i], 1-MARKS[i]); end; +function PopulateRightMarks(AOwner: TComponent): TListChartSource; +const + SIGMA = 1; + MEAN = 0; +var + i: Integer; + z: Double; +begin + Result := TListChartSource.Create(AOwner); + for i := -6 to -1 do + begin + z := NormalDist(i); //*SIGMA, MEAN, SIGMA); + Result.Add(z, z, 'μ -' + IntToStr(-i) + ' σ'); + end; + Result.Add(0.5, 0.5, 'Mean (μ)'); + for i := 1 to 6 do + begin + z := NormalDist(i); //*SIGMA, MEAN, SIGMA); + Result.Add(z, z, 'μ + ' + IntToStr(i) + ' σ'); + end; +end; + + +{ TNormalityFrm } + procedure TNormalityFrm.Calc_Lilliefors(const AData: DblDyneVec; out ASkew, AKurtosis, AStat: Double; out AConclusion: String); var @@ -222,13 +189,7 @@ begin // Obtain the test statistic for i := 1 to n1 do - begin - F1 := Norm(x[i]); - if x[i] >= 0 then - fval[i] := 1.0 - (F1 / 2.0) - else - fval[i] := F1 / 2.0; - end; + fval[i] := NormalDist(x[i]); // Cumulative proportions jval[1] := freq[1] / n; @@ -294,20 +255,94 @@ begin end; +procedure TNormalityFrm.CloseBtnClick(Sender: TObject); +begin + Close; +end; + + +procedure TNormalityFrm.FormActivate(Sender: TObject); +var + w: Integer; +begin + if FAutoSized then + exit; + + w := MaxValue([ResetBtn.Width, ComputeBtn.Width, CloseBtn.Width]); + ResetBtn.Constraints.MinWidth := w; + ComputeBtn.Constraints.MinWidth := w; + CloseBtn.Constraints.MinWidth := w; + + ParamsPanel.Constraints.MinWidth := Max( + 3*w + 2*CloseBtn.BorderSpacing.Left, + Max(Label1.Width, Label2.Width) + VarInBtn.Width + 2*VarList.BorderSpacing.Right); + ParamsPanel.Constraints.MinHeight := VarOutBtn.Top + VarOutBtn.Height + + VarOutBtn.BorderSpacing.Bottom + Bevel1.Height + Panel1.Height + + Panel1.BorderSpacing.Top; + + Constraints.MinWidth := ParamsPanel.Constraints.MinWidth + 300; + Constraints.MinHeight := ParamsPanel.Constraints.MinHeight + 2*ParamsPanel.BorderSpacing.Top; + + Position := poDesigned; + FAutoSized := True; +end; + +procedure TNormalityFrm.FormCreate(Sender: TObject); +var + T : TChartAxisTransformations; + InvNormDistTransform: TAxisTransform; + L: TListChartSource; + i: Integer; +begin + Assert(OS3MainFrm <> nil); + + InitForm(self); + + FReportFrame := TReportFrame.Create(self); + FReportFrame.Parent := ReportPage; + FReportFrame.Align := alClient; + + FChartFrame := TChartFrame.Create(self); + FChartFrame.Parent := ChartPage; + FChartFrame.Align := alClient; + FChartFrame.Chart.BottomAxis.Intervals.MaxLength := 80; + FChartFrame.Chart.BottomAxis.Intervals.MinLength := 30; + + T := TChartAxisTransformations.Create(FChartFrame); + InvNormDistTransform := TCumulNormDistrAxisTransform.Create(T); + InvNormDistTransform.Transformations := T; + + with FChartFrame.Chart.LeftAxis do + begin + Transformations := T; + + Marks.Source := PopulateLeftMarks(FChartFrame.Chart); + Marks.Style := smsValue; + Marks.OverlapPolicy := opHideNeighbour; + end; + + with TChartAxis(FChartFrame.Chart.AxisList.Add) do + begin + Alignment := calRight; + Transformations := T; + Marks.Source := PopulateRightMarks(FChartFrame.Chart); + Marks.Style := smsLabel; + Marks.TextFormat := tfHTML; + Grid.Color := clSilver; + Grid.Style := psSolid; + end; + + Reset; +end; + + procedure TNormalityFrm.ComputeBtnClick(Sender: TObject); var w: Double = 0.0; pw: Double = 0.0; - skew, kurtosis: double; - mean, variance, stddev, deviation, devsqr, M2, M3, M4: double; - i, j, n, n1: integer; + skew, kurtosis, D: double; + i, j, n: integer; data: DblDyneVec = nil; - //a - z, x: DblDyneVec; - freq: IntDyneVec; - fval, jval, DP: DblDyneVec; - F1, DPP, D, A0, C1, D15, D10, D05, D025, t2: double; - init : boolean; conclusion: String; lReport: TStrings; begin @@ -316,15 +351,7 @@ begin exit; n := Length(data) - 1; // Subtract 1 because of unused element at index 0 - (* -// SetLength(a, n+1); // +1 because all arrays are considered to begin at 1 here - SetLength(freq, n+1); - SetLength(z, n+1); - SetLength(x, n+1); - SetLength(fval, n+1); - SetLength(jval, n+1); - SetLength(DP, n+1); - *) + // Sort into ascending order for i := 1 to n - 1 do for j := i + 1 to n do @@ -333,117 +360,11 @@ begin if not Calc_ShapiroWilks(data, w, pw) then exit; - (* - // Call Shapiro-Wilks function - init := false; - swilk(init, data, n, n1, n2, a, w, pw, ier); - if ier <> 0 then - begin - ErrorMsg('Error encountered: ' + IntToStr(ier)); - Cleanup; - exit; - end; - *) - - { - WEdit.Text := Format('%.4f', [w]); - ProbEdit.Text := Format('%.4f', [pw]); - } // Now do Lilliefors Calc_Lilliefors(data, skew, kurtosis, D, conclusion); - (* - // Get unique scores and their frequencies - n1 := 1; - i := 1; - freq[1] := 1; - x[1] := data[1]; - repeat - for j := i + 1 to n do - begin - if data[j] = x[n1] then freq[n1] := freq[n1] + 1; - end; - i := i + freq[n1]; - if i <= n then - begin - n1 := n1 + 1; - x[n1] := data[i]; - freq[n1] := 1; - end; - until i > n; - - // Now get skew and kurtosis of scores - mean := 0.0; - variance := 0.0; - for i := 1 to n do - begin - mean := mean + data[i]; - variance := variance + (data[i] * data[i]); - end; - variance := variance - (mean * mean) / n; - variance := variance / (n - 1); - stddev := sqrt(variance); - mean := mean / n; - - // Obtain skew, kurtosis and z scores - M2 := 0.0; - M3 := 0.0; - M4 := 0.0; - for i := 1 to n do - begin - deviation := data[i] - mean; - devsqr := deviation * deviation; - M2 := M2 + devsqr; - M3 := M3 + (deviation * devsqr); - M4 := M4 + (devsqr * devsqr); - z[i] := (data[i] - mean) / stddev; - end; - for i := 1 to n1 do x[i] := (x[i] - mean) / stddev; - skew := (n * M3) / ((n - 1) * (n - 2) * stddev * variance); - kurtosis := (n * (n + 1) * M4) - (3 * M2 * M2 * (n - 1)); - kurtosis := kurtosis /( (n - 1) * (n - 2) * (n - 3) * (variance * variance) ); - //SkewEdit.Text := Format('%.3f', [skew]); - //KurtosisEdit.Text := Format('%.3f', [kurtosis]); - - // Obtain the test statistic - for i := 1 to n1 do - begin - F1 := Norm(x[i]); - if x[i] >= 0 then - fval[i] := 1.0 - (F1 / 2.0) - else - fval[i] := F1 / 2.0; - end; - - // Cumulative proportions - jval[1] := freq[1] / n; - for i := 2 to n1 do jval[i] := jval[i-1] + freq[i] / n; - for i := 1 to n1 do DP[i] := abs(jval[i] - fval[i]); - - // Sort DP - for i := 1 to n1-1 do - for j := i+1 to n1 do - if DP[j] < DP[i] then - Exchange(DP[i], DP[j]); - - DPP := DP[n1]; - D := DPP; - //StatEdit.Text := Format('%.3f', [D]); - A0 := sqrt(n); - C1 := A0 - 0.01 + (0.85 / A0); - D15 := 0.775 / C1; - D10 := 0.819 / C1; - D05 := 0.895 / C1; - D025 := 0.995 / C1; - t2 := D; - if t2 > D025 then conclusion := 'Strong evidence against normality.'; - if ((t2 <= D025) and (t2 > D05)) then conclusion := 'Sufficient evidence against normality.'; - if ((t2 <= D05) and (t2 > D10)) then conclusion := 'Suggestive evidence against normality.'; - if ((t2 <= D10) and (t2 > D15)) then conclusion := 'Little evidence against normality.'; - if (t2 <= D15) then conclusion := 'No evidence against normality.'; - *) - + // Print results to report frame lReport := TStringList.Create; try lReport.Add('NORMALITY TESTS FOR '+ TestVarEdit.Text); @@ -463,27 +384,20 @@ begin lReport.Free; end; + // Probability plot PlotData(data); end; -function TNormalityFrm.Norm(z : double) : double; -var - p: double; -begin - z := abs(z); - p := 1.0 + z * (0.04986735 + z * (0.02114101 + z * (0.00327763 + - z * (0.0000380036 + z * (0.0000488906 + z * 0.000005383))))); - p := p * p; - p := p * p; - p := p * p; - Result := 1.0 / (p * p); -end; +{ Plots the cumulative probability of the data onto the probability grid + prepared in FormCreate } procedure TNormalityFrm.PlotData(AData: DblDyneVec); var i, n: Integer; ser: TChartSeries; begin + FChartFrame.Clear; + n := Length(AData) - 1; // take care of unused value at index 0 ser := FChartFrame.PlotXY(ptSymbols, nil, nil, nil, nil, '', clBlack); ser.AxisIndexX := 1; @@ -563,7 +477,6 @@ begin end; - procedure TNormalityFrm.UpdateBtnStates; begin VarInBtn.Enabled := (VarList.ItemIndex > -1) and (TestVarEdit.Text = '');