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

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
+  case PlanFromCombo() of
-    0: Plan1;
+    1: Plan1;
-    1: Plan2;
+    2: Plan2;
-    2: Plan3;
+    3: Plan3;
-    3: Plan4;
+    4: Plan4;
-    4: Plan5;
+    5: Plan5;
-    5: Plan6;
+    6: Plan6;
-    6: Plan7;
+    7: Plan7;
-    7: Plan9;
+    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;
+      cellcnts[i, j] := 0;
-      celltotals[i,j] := 0.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]));
+    // Data values in cols A, B, C are assumed to be 1-based indices! --> subtract 1
-    col := Round(StrTofloat(OS3MainFrm.DataGrid.Cells[Bcol, i]));
+    row := Round(StrToFloat(OS3MainFrm.DataGrid.Cells[Acol, i])) - 1;
-    slice := Round(StrToFloat(OS3MainFrm.DataGrid.Cells[Ccol, i]));
+    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;
+    cellcnts[row, col] := cellcnts[row, col] + 1;
-    celltotals[row-1,col-1] := celltotals[row-1,col-1] + data;
+    celltotals[row, col] := celltotals[row, col] + data;
-    Ctotals[slice-1] := Ctotals[slice-1] + data;
+    Ctotals[slice] := Ctotals[slice] + data;
-    sumxsqr := sumxsqr + (data * 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[i, p] := celltotals[i, p] + celltotals[i, j];
-      celltotals[p,j] := celltotals[p,j] + celltotals[i,j];
+      celltotals[p, j] := celltotals[p, j] + celltotals[i, j];
-      sumABCsqr := sumABCsqr + (celltotals[i,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
+  // Collapse slices into group x a matrix
-  // result is the group times A matrix with BC cells containing n cases each
+  // 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 [2, 3, 4, 7, 9];
-  DInBtn.Visible := APlan in [1, 2, 3, 6, 8];
   DOutBtn.Visible := DInBtn.Visible;
   DCodeLabel.Visible := DInBtn.Visible;
   DCodeEdit.Visible := DInBtn.Visible;
 
-//  GrpInBtn.Visible := APlan in [5, 6, 7, 9];
+  GrpInBtn.Visible := APlan in [5, 6, 7, 9];
-  GrpInBtn.Visible := APlan in [4, 5, 6, 7];
   GrpOutBtn.Visible := GrpInBtn.Visible;
   GrpCodeLabel.Visible := GrpInBtn.Visible;
   GrpCodeEdit.Visible := GrpInBtn.Visible;

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;