From 7800484f54517c959618346685d0097b8271bcad Mon Sep 17 00:00:00 2001 From: wp_xxyyzz Date: Sat, 5 Dec 2020 18:29:18 +0000 Subject: [PATCH] LazStats: Fix crash of Latin-Square plan #7 git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7924 8e941d3f-bd1b-0410-a28a-d453659cc2b4 --- .../analysis/comparisons/latinsqrsunit.pas | 110 +++++++++++------- .../lazstats/source/units/mathunit.pas | 5 +- 2 files changed, 69 insertions(+), 46 deletions(-) diff --git a/applications/lazstats/source/forms/analysis/comparisons/latinsqrsunit.pas b/applications/lazstats/source/forms/analysis/comparisons/latinsqrsunit.pas index e65078fa9..2255abe38 100644 --- a/applications/lazstats/source/forms/analysis/comparisons/latinsqrsunit.pas +++ b/applications/lazstats/source/forms/analysis/comparisons/latinsqrsunit.pas @@ -75,6 +75,7 @@ type // procedure Plan8; procedure Plan9; + function PlanFromCombo: Integer; procedure PrepareForPlan(APlan: Integer); procedure ResetPlan; @@ -220,15 +221,15 @@ end; procedure TLatinSqrsForm.Compute; begin - case PlanCombo.ItemIndex of - 0: Plan1; - 1: Plan2; - 2: Plan3; - 3: Plan4; - 4: Plan5; - 5: Plan6; - 6: Plan7; - 7: Plan9; + case PlanFromCombo() of + 1: Plan1; + 2: Plan2; + 3: Plan3; + 4: Plan4; + 5: Plan5; + 6: Plan6; + 7: Plan7; + 9: Plan9; end; end; @@ -334,7 +335,7 @@ end; procedure TLatinSqrsForm.PlanComboChange(Sender: TObject); begin PlanPanel.Visible := true; - PrepareForPlan(PlanCombo.ItemIndex); + PrepareForPlan(PlanFromCombo); end; @@ -388,7 +389,7 @@ begin p := rangeA; // set up an array for cell counts and for cell sums and marginal sums - SetLength(cellcnts, rangeA+1, rangeB+1); + SetLength(cellcnts, rangeA+1, rangeB+1); // wp: why +1? SetLength(celltotals, rangeA+1, rangeB+1); SetLength(Ctotals, rangeC+1); SetLength(Design, rangeA, rangeB); @@ -397,8 +398,8 @@ begin for i := 0 to rangeA do for j := 0 to rangeB do begin - cellcnts[i,j] := 0; - celltotals[i,j] := 0.0; + cellcnts[i, j] := 0; + celltotals[i, j] := 0.0; end; for i := 0 to rangeC-1 do Ctotals[i] := 0; @@ -417,21 +418,22 @@ begin term6 := 0.0; GrandMean := 0.0; - // Read in the data + // Read in the data for i := 1 to NoCases do begin - row := Round(StrToFloat(OS3MainFrm.DataGrid.Cells[Acol, i])); - col := Round(StrTofloat(OS3MainFrm.DataGrid.Cells[Bcol, i])); - slice := Round(StrToFloat(OS3MainFrm.DataGrid.Cells[Ccol, i])); + // Data values in cols A, B, C are assumed to be 1-based indices! --> subtract 1 + row := Round(StrToFloat(OS3MainFrm.DataGrid.Cells[Acol, i])) - 1; + col := Round(StrTofloat(OS3MainFrm.DataGrid.Cells[Bcol, i])) - 1; + slice := Round(StrToFloat(OS3MainFrm.DataGrid.Cells[Ccol, i])) - 1; data := StrToFloat(OS3MainFrm.DataGrid.Cells[DataCol, i]); - cellcnts[row-1,col-1] := cellcnts[row-1,col-1] + 1; - celltotals[row-1,col-1] := celltotals[row-1,col-1] + data; - Ctotals[slice-1] := Ctotals[slice-1] + data; - sumxsqr := sumxsqr + (data * data); + cellcnts[row, col] := cellcnts[row, col] + 1; + celltotals[row, col] := celltotals[row, col] + data; + Ctotals[slice] := Ctotals[slice] + data; + sumxsqr := sumxsqr + sqr(data); GrandMean := GrandMean + data; end; - // check for equal cell counts + // Check for equal cell counts for i := 0 to p-1 do for j := 0 to p-1 do begin @@ -442,14 +444,14 @@ begin end; end; - // calculate values + // Calculate values for i := 0 to p - 1 do // get row and column sums begin - for j := 0 to p-1 do + for j := 0 to p - 1 do begin - celltotals[i,p] := celltotals[i,p] + celltotals[i,j]; - celltotals[p,j] := celltotals[p,j] + celltotals[i,j]; - sumABCsqr := sumABCsqr + (celltotals[i,j] * celltotals[i,j]); + celltotals[i, p] := celltotals[i, p] + celltotals[i, j]; + celltotals[p, j] := celltotals[p, j] + celltotals[i, j]; + sumABCsqr := sumABCsqr + sqr(celltotals[i,j]); end; end; @@ -536,7 +538,7 @@ begin begin cellstring := Format(' %3d ',[i+1]); for j := 0 to p - 1 do - cellstring := cellstring + Format('%5s', [Design[i,j]]); + cellstring := cellstring + Format('%5s', [Design[i, j]]); lReport.Add(cellstring); end; cellstring := '----------'; @@ -3141,14 +3143,14 @@ begin exit; end; - // collapse slices into group x a matrix - // result is the group times A matrix with BC cells containing n cases each + // Collapse slices into group x a matrix + // Result is the group times A matrix with BC cells containing n cases each for i := 0 to p-1 do // group for j := 0 to p-1 do // factor a for k := 0 to n-1 do // factor c ABmat[i,j] := ABmat[i,j] + ABCmat[i,j,k]; - // get marginal totals for ABmat, GBmat and GCmat + // Get marginal totals for ABmat, GBmat and GCmat for i := 0 to p-1 do for j := 0 to p-1 do begin @@ -3160,7 +3162,7 @@ begin GCmat[i,p] := GCmat[i,p] + GCmat[i,j]; end; - // get grand total for ABmat, GBmat and GCmat + // Get grand total for ABmat, GBmat and GCmat for i := 0 to p-1 do begin ABmat[p,p] := ABmat[p,p] + ABmat[p,i]; @@ -3207,6 +3209,7 @@ begin lReport.Add(''); // get squared sum of subject's totals in each group + term7 := 0.0; for i := 0 to p-1 do // group term7 := term7 + (Subjtotals[i,n] * Subjtotals[i,n]); term7 := term7 / (n*p); // Sum G^2 sub k @@ -3223,29 +3226,35 @@ begin Subjtotals[i,n] := Subjtotals[i,n] + Subjtotals[i,j]; // get sum of squares for subjects within groups - for i := 0 to p-1 do term6 := term6 + Subjtotals[i,n]; + term6 := 0.0; + for i := 0 to p-1 do + term6 := term6 + Subjtotals[i,n]; SSsubwgrps := (term6 / p) - term7; - // get correction term and term for total sum of squares + // Get correction term and term for total sum of squares term1 := sqr(GrandMean) / (n * p * p); term2 := sumxsqr; - // get sum of squared groups for term4 of sum of squares for groups + // Get sum of squared groups for term4 of sum of squares for groups + term4 := 0.0; for i := 0 to p-1 do term4 := term4 + sqr(Grptotals[i]); term4 := term4 / (n * p); - // get sum of squared a's for term3 + // Get sum of squared a's for term3 + term3 := 0.0; for j := 0 to p-1 do // levels of a term3 := term3 + sqr(Atotals[j]); term3 := term3 / (n * p); - // get squared sum of b's (across groups) for term5 of sum of squares b + // Get squared sum of b's (across groups) for term5 of sum of squares b + term5 := 0.0; for j := 0 to p-1 do term5 := term5 + sqr(Btotals[j]); term5 := term5 / (n * p); - // get squared sum of c's (across groups) for term8 of SS for c + // Get squared sum of c's (across groups) for term8 of SS for c + term8 := 0.0; for j := 0 to p-1 do term8 := term8 + sqr(Ctotals[j]); term8 := term8 / (n * p); @@ -3258,7 +3267,7 @@ begin SSb := term5 - term1; SSc := term8 - term1; - // get sum of squared AB cells for term6 + // Get sum of squared AB cells for term6 term6 := 0.0; for i := 0 to p-1 do for j := 0 to p-1 do @@ -3291,12 +3300,14 @@ begin fa := MSa / MSerrwithin; fb := MSb / MSerrwithin; fc := MSc / MSerrwithin; - if dfab > 0 then fab := MSab / MSerrwithin; + if dfab > 0 then + fab := MSab / MSerrwithin; probgrps := probf(fgroups, dfgroups, dfsubwgrps); proba := probf(fa, dfa, dferrwithin); probb := probf(fb, dfb, dferrwithin); probc := probf(fc, dfc, dferrwithin); - probab := probf(fab, dfab, dferrwithin); + if dfab > 0 then + probab := probf(fab, dfab, dferrwithin); // show ANOVA table results lReport.Add('LATIN SQUARES REPEATED ANALYSIS Plan 7 (superimposed squares)'); @@ -4164,18 +4175,27 @@ begin end; end; + +function TLatinSqrsForm.PlanFromCombo: Integer; +begin + if PlanCombo.ItemIndex = PlanCombo.Items.Count - 1 then + // Correct for missing Plan 8 + Result := PlanCombo.ItemIndex + 2 + else + Result := PlanCombo.ItemIndex + 1; +end; + + procedure TLatinSqrsForm.PrepareForPlan(APlan: Integer); begin ResetPlan; -// DInBtn.Visible := APlan in [2, 3, 4, 7, 9]; - DInBtn.Visible := APlan in [1, 2, 3, 6, 8]; + DInBtn.Visible := APlan in [2, 3, 4, 7, 9]; DOutBtn.Visible := DInBtn.Visible; DCodeLabel.Visible := DInBtn.Visible; DCodeEdit.Visible := DInBtn.Visible; -// GrpInBtn.Visible := APlan in [5, 6, 7, 9]; - GrpInBtn.Visible := APlan in [4, 5, 6, 7]; + GrpInBtn.Visible := APlan in [5, 6, 7, 9]; GrpOutBtn.Visible := GrpInBtn.Visible; GrpCodeLabel.Visible := GrpInBtn.Visible; GrpCodeEdit.Visible := GrpInBtn.Visible; diff --git a/applications/lazstats/source/units/mathunit.pas b/applications/lazstats/source/units/mathunit.pas index 1d145999a..51cbd7691 100644 --- a/applications/lazstats/source/units/mathunit.pas +++ b/applications/lazstats/source/units/mathunit.pas @@ -407,9 +407,12 @@ begin Result := 1.0 else begin + Result := betai(0.5*df2, 0.5*df1, df2/(df2+df1*F)); // wp: this is the version from NumLib + { b1 := betai(0.5*df2, 0.5*df1, df2/(df2+df1*f)); b2 := betai(0.5*df1, 0.5*df2, df1/(df1+df2/f)); // wp here "/f" but in prev line "*f" ???? - Result := (b1 + (1.0-b2)) * 0.5; + Result := (b1 + (1.0-b2)) * 0.5; // wp: looks like average of two calc methos + } end; end;