pascal/lazarus-ccr
1
0
Fork 0
Code Issues Releases Activity

LazStats: Complete refactored NormalityUnit.

Browse Source 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7697 8e941d3f-bd1b-0410-a28a-d453659cc2b4
This commit is contained in:
wp_xxyyzz
2020-09-26 21:10:53 +00:00
parent 3c23258256
commit e353889d67
2 changed files with 127 additions and 213 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/lazstats/source/forms/analysis/descriptive/normalityunit.lfm
View File

@ -54,6 +54,7 @@ object NormalityFrm: TNormalityFrm
BorderSpacing.Left = 8 BorderSpacing.Left = 8
Caption = 'Close' Caption = 'Close'
ModalResult = 11 ModalResult = 11
OnClick = CloseBtnClick
TabOrder = 0 TabOrder = 0
end end
object ComputeBtn: TButton object ComputeBtn: TButton

339
applications/lazstats/source/forms/analysis/descriptive/normalityunit.pas
View File

@ -33,6 +33,7 @@ type
VarOutBtn: TBitBtn; VarOutBtn: TBitBtn;
Label1: TLabel; Label1: TLabel;
VarList: TListBox; VarList: TListBox;
procedure CloseBtnClick(Sender: TObject);
procedure ComputeBtnClick(Sender: TObject); procedure ComputeBtnClick(Sender: TObject);
procedure FormActivate(Sender: TObject); procedure FormActivate(Sender: TObject);
procedure FormCreate(Sender: TObject); procedure FormCreate(Sender: TObject);
@ -49,7 +50,6 @@ type
function Calc_ShapiroWilks(const AData: DblDyneVec; out W, Prob: Double): Boolean; function Calc_ShapiroWilks(const AData: DblDyneVec; out W, Prob: Double): Boolean;
procedure Calc_Lilliefors(const AData: DblDyneVec; procedure Calc_Lilliefors(const AData: DblDyneVec;
out ASkew, AKurtosis, AStat: Double; out AConclusion: String); out ASkew, AKurtosis, AStat: Double; out AConclusion: String);
function Norm(z : double): double;
procedure PlotData(AData: DblDyneVec); procedure PlotData(AData: DblDyneVec);
function PrepareData(const VarName: String): DblDyneVec; function PrepareData(const VarName: String): DblDyneVec;
procedure UpdateBtnStates; procedure UpdateBtnStates;
@ -67,84 +67,51 @@ implementation
{$R *.lfm} {$R *.lfm}
uses uses
Math, Math, MathUnit,
TAChartUtils, TATextElements, TACustomSeries, TATransformations, TAChartUtils, TATextElements, TACustomSeries, TATransformations,
TACustomSource, TASources, TAChartAxisUtils, TAChartAxis, TACustomSource, TASources,
Utils; Utils;
{ TNormalityFrm }
procedure TNormalityFrm.FormActivate(Sender: TObject); function PopulateLeftMarks(AOwner: TComponent): TListChartSource;
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);
const const
MARKS: array[0..23] of double = (0.001, 0.002, 0.005, 0.01, 0.02, 0.03, 0.05, 0.07, MARKS: array[0..12] of double = (
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.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);
0.995, 0.998, 0.999);
var var
T : TChartAxisTransformations;
InvNormDistTransform: TAxisTransform;
L: TListChartSource;
i: Integer; i: Integer;
begin begin
Assert(OS3MainFrm <> nil); Result := TListChartSource.Create(AOwner);
for i := 0 to High(MARKS) do Result.Add(MARKS[i], MARKS[i]);
InitForm(self); Result.Add(0.5, 0.5);
for i := High(Marks) downto 0 do Result.Add(1-MARKS[i], 1-MARKS[i]);
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;
end; 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, '&mu; -' + IntToStr(-i) + ' &sigma;');
end;
Result.Add(0.5, 0.5, 'Mean (&mu;)');
for i := 1 to 6 do
begin
z := NormalDist(i); //*SIGMA, MEAN, SIGMA);
Result.Add(z, z, '&mu; + ' + IntToStr(i) + ' &sigma;');
end;
end;
{ TNormalityFrm }
procedure TNormalityFrm.Calc_Lilliefors(const AData: DblDyneVec; procedure TNormalityFrm.Calc_Lilliefors(const AData: DblDyneVec;
out ASkew, AKurtosis, AStat: Double; out AConclusion: String); out ASkew, AKurtosis, AStat: Double; out AConclusion: String);
var var
@ -222,13 +189,7 @@ begin
// Obtain the test statistic // Obtain the test statistic
for i := 1 to n1 do for i := 1 to n1 do
begin fval[i] := NormalDist(x[i]);
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 // Cumulative proportions
jval[1] := freq[1] / n; jval[1] := freq[1] / n;
@ -294,20 +255,94 @@ begin
end; 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); procedure TNormalityFrm.ComputeBtnClick(Sender: TObject);
var var
w: Double = 0.0; w: Double = 0.0;
pw: Double = 0.0; pw: Double = 0.0;
skew, kurtosis: double; skew, kurtosis, D: double;
mean, variance, stddev, deviation, devsqr, M2, M3, M4: double; i, j, n: integer;
i, j, n, n1: integer;
data: DblDyneVec = nil; 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; conclusion: String;
lReport: TStrings; lReport: TStrings;
begin begin
@ -316,15 +351,7 @@ begin
exit; exit;
n := Length(data) - 1; // Subtract 1 because of unused element at index 0 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 // Sort into ascending order
for i := 1 to n - 1 do for i := 1 to n - 1 do
for j := i + 1 to n do for j := i + 1 to n do
@ -333,117 +360,11 @@ begin
if not Calc_ShapiroWilks(data, w, pw) then if not Calc_ShapiroWilks(data, w, pw) then
exit; 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 // Now do Lilliefors
Calc_Lilliefors(data, skew, kurtosis, D, conclusion); Calc_Lilliefors(data, skew, kurtosis, D, conclusion);
(* // Print results to report frame
// 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.';
*)
lReport := TStringList.Create; lReport := TStringList.Create;
try try
lReport.Add('NORMALITY TESTS FOR '+ TestVarEdit.Text); lReport.Add('NORMALITY TESTS FOR '+ TestVarEdit.Text);
@ -463,27 +384,20 @@ begin
lReport.Free; lReport.Free;
end; end;
// Probability plot
PlotData(data); PlotData(data);
end; 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); procedure TNormalityFrm.PlotData(AData: DblDyneVec);
var var
i, n: Integer; i, n: Integer;
ser: TChartSeries; ser: TChartSeries;
begin begin
FChartFrame.Clear;
n := Length(AData) - 1; // take care of unused value at index 0 n := Length(AData) - 1; // take care of unused value at index 0
ser := FChartFrame.PlotXY(ptSymbols, nil, nil, nil, nil, '', clBlack); ser := FChartFrame.PlotXY(ptSymbols, nil, nil, nil, nil, '', clBlack);
ser.AxisIndexX := 1; ser.AxisIndexX := 1;
@ -563,7 +477,6 @@ begin
end; end;
procedure TNormalityFrm.UpdateBtnStates; procedure TNormalityFrm.UpdateBtnStates;
begin begin
VarInBtn.Enabled := (VarList.ItemIndex > -1) and (TestVarEdit.Text = ''); VarInBtn.Enabled := (VarList.ItemIndex > -1) and (TestVarEdit.Text = '');