LazStats: Fix individuals x-bar chart to agree with JMP software. Nicer labeling of horizontal lines in the SPC charts. Calculate and display Cp, Cpk and Cpm values.

git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7665 8e941d3f-bd1b-0410-a28a-d453659cc2b4

applications/lazstats/source/forms/analysis/statistical_process_control/basicspcunit.pas

@@ -128,6 +128,8 @@ const
   TOP_MARGIN = 150;
   BOTTOM_MARGIN = 200;
 
+  FORMAT_MASK = '0.000';
+
 
 { TBasicSPCForm }
 
@@ -323,7 +325,7 @@ const
   TARGET_COLOR = clBlue;
   CL_COLOR = clRed;
   SPEC_COLOR = clGreen;
-  CL_STYLE = psDot;
+  CL_STYLE = psDash;
   SPEC_STYLE = psSolid;
 var
   ser: TChartSeries;
@@ -349,19 +351,19 @@ begin
   if not IsNaN(GrandMean) then
   begin
     FChartFrame.HorLine(GrandMean, clRed, psSolid, AGrandMeanTitle);
-    rightLabels.Add(GrandMean, GrandMean, AGrandMeanTitle);
+    rightLabels.Add(GrandMean, GrandMean, AGrandMeanTitle + '=' + FormatFloat(FORMAT_MASK, GrandMean));
   end;
 
   if not IsNaN(UCL) then
   begin
     FChartFrame.HorLine(UCL, CL_COLOR, CL_STYLE, 'UCL/LCL');
-    rightLabels.Add(UCL, UCL, 'UCL');
+    rightLabels.Add(UCL, UCL, 'UCL=' + FormatFloat(FORMAT_MASK, UCL));
   end;
 
   if not IsNaN(LCL) then
   begin
     FChartFrame.HorLine(LCL, CL_COLOR, CL_STYLE, '');
-    rightLabels.Add(UCL, LCL, 'LCL');
+    rightLabels.Add(UCL, LCL, 'LCL=' + FormatFloat(FORMAT_MASK, LCL));
   end;
 
   if not IsNan(UpperSpec) then
@@ -371,12 +373,12 @@ begin
     else
       s := 'Upper/Lower Spec';
     FChartFrame.HorLine(UpperSpec, SPEC_COLOR, SPEC_STYLE, s);
-    rightLabels.Add(UpperSpec, UpperSpec, 'Upper Spec');
+    rightLabels.Add(UpperSpec, UpperSpec, 'USL=' + FormatFloat(FORMAT_MASK, UpperSpec));
   end;
 
   if not IsNaN(TargetSpec) then begin
     FChartFrame.HorLine(TargetSpec, TARGET_COLOR, psSolid, 'Target');
-    rightLabels.Add(TargetSpec, TargetSpec, 'Target');
+    rightLabels.Add(TargetSpec, TargetSpec, 'Target=' + FormatFloat(FORMAT_MASK, TargetSpec));
   end;
 
   if not IsNaN(LowerSpec) then
@@ -387,8 +389,10 @@ begin
       s := 'Upper/Lower Spec';
     constLine := FChartFrame.HorLine(LowerSpec, SPEC_COLOR, SPEC_STYLE, s);
     constLine.Legend.Visible := IsNaN(UpperSpec);
-    rightLabels.Add(LowerSpec, LowerSpec, 'Lower Spec');
+    rightLabels.Add(LowerSpec, LowerSpec, 'LSL=' + FormatFloat(FORMAT_MASK, LowerSpec));
   end;
+
+  FChartFrame.Chart.Legend.Visible := false;
 end;
 
 

applications/lazstats/source/forms/analysis/statistical_process_control/xbarchartunit.lfm

@@ -1,19 +1,19 @@
 inherited XBarChartForm: TXBarChartForm
   Left = 572
-  Height = 431
+  Height = 476
   Top = 215
   HelpType = htKeyword
   HelpKeyword = 'html/XBarChart.htm'
   Caption = 'X-Bar Control Chart'
-  ClientHeight = 431
+  ClientHeight = 476
   OnActivate = FormActivate
   inherited SpecsPanel: TPanel
-    Height = 431
+    Height = 476
     Width = 379
-    ClientHeight = 431
+    ClientHeight = 476
     ClientWidth = 379
     inherited ButtonPanel: TPanel
-      Top = 389
+      Top = 434
       Width = 379
       ClientWidth = 379
       TabOrder = 5
@@ -34,7 +34,7 @@ inherited XBarChartForm: TXBarChartForm
       end
     end
     inherited VarList: TListBox
-      Height = 356
+      Height = 401
       Width = 169
     end
     inherited GroupLabel: TLabel
@@ -121,7 +121,7 @@ inherited XBarChartForm: TXBarChartForm
       AnchorSideRight.Control = MeasEdit
       AnchorSideRight.Side = asrBottom
       Left = 177
-      Height = 103
+      Height = 128
       Top = 293
       Width = 202
       Anchors = [akTop, akLeft, akRight]
@@ -129,7 +129,7 @@ inherited XBarChartForm: TXBarChartForm
       BorderSpacing.Top = 12
       BorderSpacing.Bottom = 8
       Caption = 'Show...'
-      ClientHeight = 83
+      ClientHeight = 108
       ClientWidth = 198
       TabOrder = 4
       object UpperSpecChk: TCheckBox
@@ -218,26 +218,38 @@ inherited XBarChartForm: TXBarChartForm
         TabOrder = 5
         Text = 'TargetSpecEdit'
       end
+      object ZonesChk: TCheckBox
+        AnchorSideLeft.Control = TargetChk
+        AnchorSideTop.Control = TargetSpecEdit
+        AnchorSideTop.Side = asrBottom
+        Left = 12
+        Height = 19
+        Top = 81
+        Width = 52
+        BorderSpacing.Bottom = 8
+        Caption = 'Zones'
+        TabOrder = 6
+      end
     end
   end
   inherited SpecsSplitter: TSplitter
     Left = 390
-    Height = 431
+    Height = 476
   end
   inherited PageControl: TPageControl
     Left = 398
-    Height = 419
+    Height = 464
     Width = 523
     inherited ReportPage: TTabSheet
-      ClientHeight = 391
+      ClientHeight = 436
       ClientWidth = 515
       inherited Panel1: TPanel
-        Height = 351
+        Height = 396
         Width = 503
-        ClientHeight = 347
+        ClientHeight = 392
         ClientWidth = 499
         inherited ReportMemo: TMemo
-          Height = 339
+          Height = 384
           Width = 491
         end
       end

applications/lazstats/source/forms/analysis/statistical_process_control/xbarchartunit.pas

@@ -1,3 +1,20 @@
+{ This unit was checked against the commercial statistical package JMP and
+  creates correct results.
+
+  Data file for testing: "boltsize.laz"
+  Group variable:         LotNo
+  Selected variable:      BoltLngth
+
+  The original LazStats help files suggest
+    Upper Spec Level      20.05
+    Lower Spec Level      19.95
+    Target Spec           20.00
+  but this would indicate a very poor process. Better values:
+    Upper Spec Level      21.00
+    Lower Spec Level      19.00
+    Target Spec           20.00
+}
+
 unit XBarChartUnit;
 
 {$mode objfpc}{$H+}
@@ -6,13 +23,15 @@ interface
 
 uses
   Classes, SysUtils, Forms, Controls, Graphics, Dialogs, ComCtrls,
-  ExtCtrls, StdCtrls, BasicSPCUnit;
+  ExtCtrls, StdCtrls, Buttons, PrintersDlgs,
+  Globals, BasicSPCUnit;
 
 type
 
   { TXBarChartForm }
 
   TXBarChartForm = class(TBasicSPCForm)
+    ZonesChk: TCheckBox;
     LevelOptns: TGroupBox;
     LowerSpecChk: TCheckBox;
     LowerSpecEdit: TEdit;
@@ -23,8 +42,13 @@ type
     UpperSpecEdit: TEdit;
     XSigmaEdit: TEdit;
     procedure FormActivate(Sender: TObject);
+  private
+    FAveStdDev: Double;
   protected
     procedure Compute; override;
+    procedure PlotMeans(ATitle, AXTitle, AYTitle, ADataTitle, AGrandMeanTitle: String;
+      const Groups: StrDyneVec; const Means: DblDyneVec;
+      UCL, LCL, GrandMean, TargetSpec, LowerSpec, UpperSpec: double); override;
     procedure Reset; override;
     function Validate(out AMsg: String; out AControl: TWinControl): Boolean; override;
   end;
@@ -37,7 +61,7 @@ implementation
 
 uses
   Math,
-  Globals, Utils, MainUnit, DataProcs;
+  Utils, MainUnit, DataProcs;
 
 {$R *.lfm}
 
@@ -80,15 +104,19 @@ var
   upperSpec: Double = NaN;
   lowerSpec: Double = NaN;
   targetSpec: Double = NaN;
+  Cp: Double = NaN;
+  Cpk: Double = NaN;
+  Cpm: Double = NaN;
   ColNoSelected: IntDyneVec = nil;
   groups: StrDyneVec = nil;
   means: DblDyneVec = nil;
   stdDev: DblDyneVec = nil;
   count: IntDyneVec = nil;
-  numGrps, grpIndex, totalNumCases, grpSize: Integer;
+  numValues, numGrps, grpIndex, grpSize: Integer;
   grp: String;
-  X, Xsq: Double;
+  X, Xsq, prevX: Double;
-  sigma, aveStdDev, UCL, LCL, grandMean, grandSD, SEMean, C4Value: Double;
+  sigma, UCL, LCL, grandMean, grandSD, SEMean: Double;
+  individualsChart: Boolean;
   lReport: TStrings;
 begin
   if GroupEdit.Text <> '' then
@@ -96,10 +124,12 @@ begin
     SetLength(ColNoSelected, 2);
     ColNoSelected[0] := GrpVar;
     ColNoSelected[1] := MeasVar;
+    individualsChart := false;
   end else
   begin
     SetLength(ColNoSelected, 1);
     ColNoSelected[0] := MeasVar;
+    individualsChart := true;
   end;
 
   if UpperSpecChk.Checked and (UpperSpecEdit.Text <> '') then
@@ -117,71 +147,83 @@ begin
     else raise Exception.Create('Sigma case not handled.');
   end;
 
-  if GroupEdit.Text = '' then
+  if individualsChart then
     SetLength(groups, NoCases)
   else
     groups := GetGroups;
   numGrps := Length(groups);
 
   SetLength(means, numGrps);
-  SetLength(count, numGrps);
   SetLength(stddev, numGrps);
-  SEMean := 0.0;
   grandMean := 0.0;
-  totalNumCases := 0;
+  grandSD := 0.0;
+  numValues := 0;
 
-  // calculate group means, grand mean, group sd's, seMean
+  // calculate group means, grand mean, group std devs, seMean
+  if IndividualsChart then
+  begin
+    // x-bar chart of individual measurements, no groups
+    SetLength(count, 0);    // not needed, count is always 1
+    prevX := NaN;
     for i := 1 to NoCases do
     begin
       if not GoodRecord(i, Length(ColNoSelected), ColNoSelected) then continue;
-    if GroupEdit.Text = '' then
+      grpIndex := numValues;  // is counted up in this loop
-    begin
+      X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[MeasVar, i]));
-      // individuals x-bar chart
+      Xsq := X*X;
-      grpIndex := totalNumCases;
       groups[grpIndex] := IntToStr(i);
+      means[grpIndex] := means[grpIndex] + X;
+      if not IsNaN(prevX) then
+        stddev[grpIndex-1] := abs(X - prevX);  // assume std dev to be moving range;
+        // -1 --> skip empty 1st value
+      grandMean := grandMean + X;
+      grandSD := grandSD + Xsq;
+      inc(numValues);
+      prevX := X;
+    end;
   end else
   begin
     // grouped x-bar chart
+    SetLength(count, numGrps);
+    for i := 1 to NoCases do
+    begin
+      if not GoodRecord(i, Length(ColNoSelected), ColNoSelected) then continue;
       grp := Trim(OS3MainFrm.DataGrid.Cells[GrpVar, i]);
       grpIndex := IndexOfString(groups, grp);
-    end;
       X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[MeasVar, i]));
       Xsq := X*X;
       inc(count[grpIndex]);
       means[grpIndex] := means[grpIndex] + X;
       stddev[grpIndex] := stddev[grpIndex] + Xsq;
       grandMean := grandMean + X;
-    SEMean := SEMean + Xsq;
+      grandSD := grandSD + Xsq;
-    inc(totalNumCases);
+      inc(numValues);
+    end;
   end;
 
-  SEMean := SEMean - sqr(grandMean) / totalNumCases;
+  grandSD := grandSD - sqr(grandMean) / numValues;
-  SEMean := sqrt(SEMean / (totalNumCases - 1));
+  grandSD := sqrt(grandSD / (numValues - 1));
-  grandSD := SEMean;
+  SEMean := grandSD / sqrt(numValues);
-  SEMean := SEMean / sqrt(totalNumCases);
+  grandMean := grandMean / numValues;
-  grandMean := grandMean / totalNumCases;
 
-  if (GroupEdit.Text = '') then
+  if individualsChart then
   begin
     // Individuals chart
     grpSize := 1;
-    SetLength(means, totalNumCases);
+    SetLength(means, numValues);
-    SetLength(stddev, totalNumCases);
+    Setlength(count, numValues);
-    Setlength(count, totalNumCases);
+    SetLength(stddev, numValues-1);  // -1 for the missing 1st value
-    for i := 0 to totalNumCases-1 do
+    FAveStdDev := 0;
-      stddev[i] := SEMean;
+    for i := 0 to High(stddev) do
-    aveStdDev := NaN;
+      FAveStdDev := FAveStdDev + stdDev[i];
-    if totalNumCases <= 25 then
+    FAveStdDev := FAveStdDev / Length(stddev) / 1.128;    // 1.128 is the value of d2 fo n = 2.
-      C4Value := C4[totalNumCases]
+    UCL := grandMean + sigma * FAveStdDev;
-    else
+    LCL := grandMean - sigma * FAveStdDev;
-      C4Value := 1.0;
-    UCL := grandMean + sigma * grandSD / C4Value;
-    LCL := grandMean - sigma * grandSD / C4Value;
   end else
   begin
     // Grouped chart
 
-    // Check group size - it is assumed that all groups are equally sized
+    // Check group size first; it is assumed that all groups are equally sized
     grpSize := count[0];
     for i := 1 to numGrps-1 do
       if count[i] <> grpSize then
@@ -190,9 +232,10 @@ begin
         exit;
       end;
 
-    aveStdDev := aveStdDev + stdDev[i];
+    SetLength(means, numGrps);
-
+    Setlength(count, numGrps);
-    aveStdDev := 0;
+    SetLength(stddev, numGrps);
+    FAveStdDev := 0;
     for i := 0 to numGrps-1 do
     begin
       if count[i] = 0 then
@@ -207,37 +250,68 @@ begin
         begin
           stddev[i] := stddev[i] - sqr(means[i]) / count[i];
           stddev[i] := stddev[i] / (count[i] - 1);   // Variance of group i
-          aveStdDev := aveStdDev + stdDev[i];       // Sum of variances
+          FAveStdDev := FAveStdDev + stdDev[i];      // Sum of variances
           stddev[i] := sqrt(stddev[i]);              // StdDev of group i
         end;
         means[i] := means[i] / count[i];
       end;
     end;
-    aveStdDev := sqrt(aveStdDev / (numGrps * grpSize));
+    FAveStdDev := sqrt(FAveStdDev / (numGrps * grpSize));
-    UCL := grandMean + sigma * aveStdDev;
+    UCL := grandMean + sigma * FAveStdDev;
-    LCL := grandMean - sigma * aveStdDev;
+    LCL := grandMean - sigma * FAveStdDev;
+
+    //UCL := grandMean + sigma * grandSD / sqrt(grpSize);  // this works, too, a bit more off of JMP than the above...
+    //LCL := grandMean - sigma * grandSD / sqrt(grpSize);
+
 //    UCL := grandMean + sigma * SEMean;    // old LazStats calculation -- does not agree with JMP software.
 //    LCL := grandMean - sigma * SEMean;
   end;
 
+  if not IsNaN(upperSpec) and not IsNaN(lowerSpec) then
+  begin
+    Cp := (upperSpec - lowerSpec) / (6* FAveStdDev);
+    Cpk := Min(UpperSpec - grandMean, grandMean - LowerSpec) / (3 * FAveStdDev);
+    if not IsNaN(targetSpec) then
+      Cpm := (upperSpec - lowerSpec) / (6 * sqrt(sqr(FAveStdDev) + sqr(grandMean - targetSpec)));
+  end;
+
     // Print results
   lReport := TStringList.Create;
   try
     lReport.Add('X BAR CHART RESULTS');
     lReport.Add('');
-    lReport.Add('Number of values:       %8d',   [totalNumCases]);
+    lReport.Add('Number of values:       %8d',   [numValues]);
     lReport.Add('Number of groups:       %8d',   [numGrps]);
     lReport.Add('Group size:             %8d',   [grpSize]);
     lReport.Add('');
     lReport.Add('Grand Mean:             %8.3f', [grandMean]);
     lReport.Add('Standard Deviation:     %8.3f', [grandSD]);
     lReport.Add('Standard Error of Mean: %8.3f', [SEMean]);
-    lReport.Add('Average Std Deviation:  %8.3f', [aveStdDev]);
+    lReport.Add('Average Std Deviation:  %8.3f', [FAveStdDev]);
-    lReport.Add('Lower Control Limit:    %8.3f', [LCL]);
     lReport.Add('Upper Control Limit:    %8.3f', [UCL]);
+    lReport.Add('Lower Control Limit:    %8.3f', [LCL]);
+    lReport.Add('');
+    if not IsNaN(targetSpec) then
+      lReport.Add('Target:                 %8.3f', [targetSpec]);
+    if not IsNaN(upperSpec) then
+      lReport.Add('Upper Spec Limit:       %8.3f', [upperSpec]);
+    if not IsNaN(lowerSpec) then
+      lReport.Add('Lower Spec Limit:       %8.3f', [lowerSpec]);
+    if not IsNaN(Cp) then
+      lReport.Add('Cp:                     %8.3f', [cp]);
+    if not IsNaN(Cpk) then
+      lReport.Add('Cpk:                    %8.3f', [Cpk]);
+    if not IsNaN(Cpm) then
+      lReport.Add('Cpm:                    %8.3f', [Cpm]);
     lReport.Add('');
     lReport.Add(' Group   Size    Mean    Std.Dev.');
     lReport.Add('-------  ----  --------  --------');
+    if individualsChart then
+    begin
+      lReport.Add  ('%7s  %4d  %8.2f', [groups[i], count[i], means[i]]);
+      for i := 1 to numGrps-1 do
+        lReport.Add('%7s  %4d  %8.2f  %8.2f', [groups[i], count[i], means[i], stddev[i-1]]);
+    end else
       for i := 0 to numGrps-1 do
         lReport.Add('%7s  %4d  %8.2f  %8.2f', [groups[i], count[i], means[i], stddev[i]]);
 
@@ -249,13 +323,42 @@ begin
   // Show graph
   PlotMeans(
     Format('x&#772; chart for "%s"', [GetFileName]),
-    GroupEdit.Text, MeasEdit.Text, 'Group means', 'Grand mean',
+    GroupEdit.Text, MeasEdit.Text, '', 'Avg',
     groups, means,
     UCL, LCL, grandmean,
     targetSpec, lowerSpec, upperSpec
   );
 end;
 
+procedure TXBarChartForm.PlotMeans(
+  ATitle, AXTitle, AYTitle, ADataTitle, AGrandMeanTitle: String;
+  const Groups: StrDyneVec; const Means: DblDyneVec;
+  UCL, LCL, GrandMean, TargetSpec, LowerSpec, UpperSpec: double);
+const
+  EPS = 1E-6;
+var
+  y: Double;
+begin
+  inherited;
+  if not ZonesChk.Checked then
+    exit;
+
+  y := GrandMean + FAveStdDev;
+  while y < UCL - EPS do
+  begin
+    FChartFrame.HorLine(y, clRed, psDot, '');
+    y := y + FAveStdDev;
+  end;
+
+  y := GrandMean - FAveStdDev;
+  while y > LCL + EPS do
+  begin
+    FChartFrame.HorLine(y, clRed, psDot, '');
+    y := y - FAveStdDev;
+  end;
+end;
+
+
 procedure TXBarChartForm.Reset;
 begin
   inherited;
@@ -266,6 +369,7 @@ begin
   UpperSpecChk.Checked := false;
   LowerSpecChk.Checked := false;
   TargetChk.Checked := false;
+  ZonesChk.Checked := false;
 end;
 
 function TXBarChartForm.Validate(out AMsg: String; out AControl: TWinControl): Boolean;