lazarus-ccr/applications/lazstats/source/forms/analysis/measurement_programs/raschunit.pas

// File for testing: itemdata2.laz
// Select the variables VAR1...VAR5

unit RaschUnit;

{$mode objfpc}{$H+}
{$WARN 6058 off : Call to subroutine "$1" marked as inline is not inlined}

interface

uses
  Classes, SysUtils, Forms, Controls, Graphics, Dialogs, StdCtrls, Buttons,
  ExtCtrls, ComCtrls, MainUnit, FunctionsLib, Globals,
  ChartFrameUnit, BasicStatsReportAndChartFormUnit;

type

  { TRaschForm }

  TRaschForm = class(TBasicStatsReportAndChartForm)
    ProxChk: TCheckBox;
    PlotItemDiffChk: TCheckBox;
    PlotScrsChk: TCheckBox;
    ItemFuncsChk: TCheckBox;
    ScoresPage: TTabSheet;
    ItemFuncsPage: TTabSheet;
    TestInfoPage: TTabSheet;
    TestInfoChk: TCheckBox;
    OptionsGroup: TGroupBox;
    InBtn: TBitBtn;
    Label2: TLabel;
    ItemList: TListBox;
    OutBtn: TBitBtn;
    Label1: TLabel;
    VarList: TListBox;
    procedure InBtnClick(Sender: TObject);
    procedure OutBtnClick(Sender: TObject);
    procedure VarListSelectionChange(Sender: TObject; {%H-}User: boolean);

  private
    procedure Analyze(const itemfail, grpfail: IntDyneVec;
      const f: IntDyneMat; var T: integer; const grppass, itempass: IntDyneVec;
      r, C1: integer; out min, max: double; const p2: DblDyneVec);
    procedure Expand(v1, v2: double; out xExpand, yExpand: Double);
    procedure FinishIt(r: integer; const i5: IntDyneVec;
      const rptbis, rbis, slope, mean: DblDyneVec;
      const itemfail: IntDyneVec; const P: DblDyneVec; AReport: TStrings);
    procedure Frequencies(C1, r: integer; const f: IntDyneMat;
      const rowtot, i5, s5: IntDyneVec; T: integer; const S: IntDyneVec;
      AReport: TStrings);
    procedure GetLogs(const L, L1, L2, g, g2, f2: DblDyneVec;
      const rowtot: IntDyneVec; k: integer; const s5, S: IntDyneVec;
      T, r, C1: integer; var v1, v2: double; AReport: TStrings);
    procedure GetScores(NoSelected: integer; const Selected: IntDyneVec;
      NoCases: integer; const f: IntDyneMat; const Mean, XSqr, SumXY: DblDyneVec;
      const S, X: IntDyneVec; out sumx, sumx2: double; var N: integer;
      AReport: TStrings);
    procedure Maxability(const expdcnt, d2, e2: DblDyneVec; const p1: DblDyneMat;
      const p2, P: DblDyneVec; C1, r: integer; const D: DblDyneMat;
      const s5: IntDyneVec; noloops: integer; AReport: TStrings);
    function MaxItem(const R1, d1: DblDyneVec; const p1, D: DblDyneMat;
      const e1, p2, P: DblDyneVec; const S, rowtot: IntDyneVec;
      T, r, C1 : integer) : double;
    procedure MaxOut(r, C1: integer; const i5, s5: IntDyneVec;
      const P, p2: DblDyneVec; AReport: TStrings);
    procedure Prox(const P,  p2: DblDyneVec; k, r, C1: integer;
      const L1: DblDyneVec; yexpand, xexpand: double; const g: DblDyneVec;
      T: integer; const rowtot, i5, s5: IntDyneVec; AReport: TStrings);
    function Reduce(k: integer; out r, T, C1: integer;
      const i5, rowtot, s5: IntDyneVec; const f: IntDyneMat; const S: IntDyneVec;
      AReport: TStrings): integer;
    procedure Slopes(const rptbis, rbis, slope: DblDyneVec; N: integer;
      sumx, sumx2: double; const sumxy: DblDyneVec; r: integer;
      const xsqr, mean: DblDyneVec);
    procedure TestFit(r, C1: integer; const f: IntDyneMat; const S: IntDyneVec;
      const P, P2: DblDyneVec; T: integer; AReport: TStrings);

    procedure PlotInfo(r: integer; const Info, A: DblDyneMat;
      const Slope, P: DblDyneVec);
    procedure Plot(const xyArray: DblDyneMat; ArraySize: integer);
    procedure PlotItems(r: integer; const i5: IntDyneVec; const P: DblDyneVec);
    procedure PlotScrs(C1: integer; const s5: IntDyneVec; const p2: DblDyneVec);
    procedure PlotTest(const TestInfo: DblDyneMat; ArraySize: integer;
      const Title: string{; Vdivisions, Hdivisions: integer});

  private
    FChartCombobox: TCombobox;
    FItemDiffsChartFrame: TChartFrame;
    FScoresChartFrame: TChartFrame;
    FItemFuncsChartFrame: TChartFrame;
    FTestInfoChartFrame: TChartFrame;
    procedure SelectItemFuncsPlot(Sender: TObject);
	
  protected
    procedure AdjustConstraints; override;
    procedure Compute; override;
    procedure UpdateBtnStates; override;
    function Validate(out AMsg: String; out AControl: TWinControl): Boolean; override;

  public
    constructor Create(AOwner: TComponent); override;
    procedure Reset; override;

  end; 

var
  RaschForm: TRaschForm;

implementation

{$R *.lfm}

uses
  Math,
  TAChartUtils, TACustomSeries,
  Utils, GridProcs;

{ TRaschForm }

constructor TRaschForm.Create(AOwner: TComponent);
begin
  inherited;

  FItemDiffsChartFrame := FChartFrame;  // already created by ancestor

  FScoresChartFrame := TChartFrame.Create(self);
  FScoresChartFrame.Parent := ScoresPage;
  FScoresChartFrame.Align := alClient;

  FItemFuncsChartFrame := TChartFrame.Create(self);
  FItemFuncsChartFrame.Parent := ItemFuncsPage;
  FItemFuncsChartFrame.Align := alClient;
  AddComboboxToToolbar(FItemFuncsChartFrame.ChartToolbar, 'Plots:', FChartCombobox);
  FChartCombobox.OnSelect := @SelectItemFuncsPlot;

  FTestInfoChartFrame := TChartFrame.Create(self);
  FTestInfoChartFrame.Parent := TestInfoPage;
  FTestInfoChartFrame.Align := alClient;
end;


procedure TRaschForm.AdjustConstraints;
begin
  inherited;
  ParamsPanel.Constraints.MinWidth := Max(
    4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left,
    OptionsGroup.Width
  );
  ParamsPanel.Constraints.MinHeight := OutBtn.Top + OutBtn.Height +
    VarList.BorderSpacing.Bottom + ButtonBevel.Height +
    CloseBtn.Borderspacing.Top + CloseBtn.Height;
end;


procedure TRaschForm.Analyze(const itemfail, grpfail: IntDyneVec;
  const f: IntDyneMat; var T: integer; const grppass, itempass: IntDyneVec;
  r, C1: integer; out min, max: double; const p2: DblDyneVec);
var
  i, j: integer;
begin
  for i := 0 to r-1 do itemfail[i] := 0;
  for j := 0 to C1-1 do grpfail[j] := 0;
  for i := 0 to r-1 do
  begin
    for j := 0 to C1-1 do
    begin
      grpfail[j] := grpfail[j] + f[i,j];
      itemfail[i] := itemfail[i] + f[i,j];
    end;
  end;
  T := 0;
  for j := 0 to C1-1 do T := T + grpfail[j];
  for j := 0 to C1-1 do grppass[j] := T - grpfail[j];
  for i := 0 to r-1 do itempass[i] := T - itemfail[i];
  min := p2[0];
  max := p2[0];
  for i := 0 to C1-1 do
  begin
    if (p2[i] < min) then min := p2[i];
    if (p2[i] > max) then max := p2[i];
  end;
end;


procedure TRaschForm.Compute;
var
   i, k1, N, C1, r, T,noloops : integer;
   sumx, sumx2, v1, v2, xexpand, yexpand, d9, min, max : double;
   X: IntDyneVec = nil;
   rowtot: IntDyneVec = nil;
   itemFail: IntDyneVec = nil;
   itemPass: IntDyneVec = nil;
   grpFail: IntDyneVec = nil;
   grpPass: IntDyneVec = nil;
   S: IntDyneVec = nil;
   s5: IntDyneVec = nil;
   i5: IntDyneVec = nil;
   f: IntDyneMat = nil;
   mean: DblDyneVec = nil;
   xsqr: DblDyneVec = nil;
   sumXY: DblDyneVec = nil;
   L: DblDyneVec = nil;
   L1: DblDyneVec = nil;
   L2: DblDyneVec = nil;
   g: DblDyneVec = nil;
   g2: DblDyneVec = nil;
   f2: DblDyneVec = nil;
   P: DblDyneVec = nil;
   p2: DblDyneVec = nil;
   R1: DblDyneVec = nil;
   d1: DblDyneVec = nil;
   d2: DblDyneVec = nil;
   e1: DblDyneVec = nil;
   e2: DblDyneVec = nil;
   expdcnt: DblDyneVec = nil;
   rptbis: DblDyneVec = nil;
   rbis: DblDyneVec = nil;
   slope: DblDyneVec = nil;
   p1: DblDyneMat = nil;
   D: DblDyneMat = nil;
   info: DblDyneMat = nil;
   A: DblDyneMat = nil;
   NoSelected : integer;
   ColNoSelected : IntDyneVec = nil;
   finished : boolean;
   error : integer;
   lReport: TStrings;
begin
   SetLength(ColNoSelected,NoVariables);
   SetLength(mean,NoVariables);
   SetLength(xsqr,NoVariables);
   SetLength(sumxy,NoVariables);
   SetLength(L,NoVariables);
   SetLength(L1,NoVariables);
   SetLength(L2,NoVariables);
   SetLength(g,NoVariables);
   SetLength(g2,NoVariables);
   SetLength(f2,NoVariables);
   SetLength(P,NoVariables);
   SetLength(p2,NoVariables);
   SetLength(R1,NoVariables);
   SetLength(d1,NoVariables);
   SetLength(e1,NoVariables);
   SetLength(expdcnt,NoVariables);
   SetLength(d2,NoVariables);
   SetLength(e2,NoVariables);
   SetLength(rptbis,NoVariables);
   SetLength(rbis,NoVariables);
   SetLength(slope,NoVariables);
   SetLength(p1,NoVariables,NoVariables);
   SetLength(D,NoVariables,NoVariables);
   SetLength(info,52,52);
   SetLength(A,52,2);
   SetLength(X,NoVariables);
   SetLength(rowtot,NoVariables);
   SetLength(itemfail,NoVariables);
   SetLength(itempass,NoVariables);
   SetLength(grpfail,NoVariables);
   SetLength(grppass,NoVariables);
   SetLength(S,NoVariables+2);
   SetLength(s5,NoVariables);
   SetLength(i5,NoVariables);
   SetLength(f,NoVariables+2,NoVariables+2);

   // Get selected variables
   NoSelected := ItemList.Items.Count;
   for i := 0 to NoSelected-1 do
     colNoSelected[i] := GetVariableIndex(OS3MainFrm.DataGrid, ItemList.Items[i]);

   //begin ( main program )
   finished := false;
   N := NoCases;
   k1 := NoSelected;

   lReport := TStringList.Create;
   try
     GetScores(NoSelected, ColNoSelected, NoCases, f, mean, xsqr, sumxy, S, X, sumx, sumx2, N, lReport);
     error := Reduce(k1, r, T, C1, i5, rowtot, s5, f, S, lReport);
     if error = 1 then
       exit;

     Frequencies(C1, r, f, rowtot, i5, s5, T, S, lReport);
     v1 := 0.0;
     v2 := 0.0;
     GetLogs(L, L1, L2, g, g2, f2, rowtot, k1, s5, S, T, r, C1, v1, v2, lReport);
     Expand(v1, v2, xexpand, yexpand);
     Prox(P, p2, k1, r, C1, L1, yexpand, xexpand, g, T, rowtot, i5, s5, lReport);
     // start iterations for the maximum-likelihood (SetLengthton-Rhapson procedure)
     // estimates
     noloops := 0;

     while (not finished) do
     begin
       d9 := MaxItem(R1, d1, p1, D, e1, p2, P, S, rowtot, T, r, C1);
       if (d9 < 0.01) then
         finished := true
       else
         Maxability(expdcnt, d2, e2, p1, p2, P, C1, r, D, s5, noloops, lReport);
           noloops := noloops + 1;
           if (noloops > 25) then
           begin
             MessageDlg('Maximum Likelihood failed to converge after 25 iterations', mtInformation, [mbOK], 0);
             finished := true;
           end;
     end;
     MaxOut(r, C1, i5, s5, P, p2, lReport);
     TestFit(r, C1, f, S, P, p2, T, lReport);
     Slopes(rptbis, rbis, slope, N, sumx, sumx2, sumxy, r, xsqr, mean);
     Analyze(itemfail, grpfail, f, T, grppass, itempass, r, C1, min, max, p2);

     if PlotItemDiffChk.Checked then PlotItems(r, i5, P);
     if PlotScrsChk.Checked then PlotScrs(C1, s5, p2);
     PlotInfo(r, info, A, slope, P);

     ChartPage.TabVisible := PlotItemDiffChk.Checked;
     ScoresPage.TabVisible := PlotScrsChk.Checked;
     ItemFuncsPage.TabVisible := ItemFuncsChk.Checked;
     TestInfopage.TabVisible := TestInfoChk.Checked;

     FinishIt(r, i5, rptbis, rbis, slope, mean, itemfail, P, lReport);

     FReportFrame.DisplayReport(lReport);

  finally
    lReport.Free;
  end;
end;


procedure TRaschForm.Expand(v1, v2: double; out xExpand, yExpand: Double);
begin
  yExpand := sqrt( (1.0 + (v2 / 2.89)) / (1.0 - (v1 * v2 / 8.35)) );
  xExpand := sqrt( (1.0 + (v1 / 2.89)) / (1.0 - (v1 * v2 / 8.35)) );
end;


procedure TRaschForm.FinishIt(r: integer; const i5: IntDyneVec;
  const rptbis, rbis, slope, mean: DblDyneVec; const itemfail: IntDyneVec;
  const P: DblDyneVec; AReport: TStrings);
var
  i: integer;
begin
  AReport.Add('');
  AReport.Add('Item Data Summary');
  AReport.Add('');
  AReport.Add('ITEM  PT.BIS.R.  BIS.R.  SLOPE   PASSED  FAILED  RASCH DIFF');
  AReport.Add('----  ---------  ------  -----   ------  ------  ----------');
             //xxx     xxxxxx   xxxxxx  xxxxx   xxxxxx   xxxx     xxxxxx
  for i := 0 to r-1 do
    AReport.Add('%3d     %6.3f   %6.3f  %5.2f   %6.2f   %4d     %6.3f', [
      i5[i], rptbis[i], rbis[i], slope[i], mean[i], itemfail[i], P[i]
    ]);
  AReport.Add('');
end;


procedure TRaschForm.Frequencies(C1, r: integer;
  const f: IntDyneMat; const rowtot, i5, s5: IntDyneVec; T: integer;
  const S: IntDyneVec; AReport: TStrings);
var
  i, j, c2, c3: integer;
  done: boolean;
  outline: string;
begin
  done := false;
  c3 := C1;
  c2 := 1;
  if (c3 > 16) then
    c3 := 16;

  while not done do
  begin
    AReport.Add('Matrix of Item Failures in Score Groups');
    outline := '   Score Group';
    for j := c2 to c3 do
        outline := outline + Format('%4d', [s5[j-1]]);
    outline := outline + '     Total';
    AReport.Add(outline);

    AReport.Add('ITEM' );
    AReport.Add('');
    for i := 1 to r do
    begin
      outline := Format('%4d          ', [i5[i-1]]);
      for j := c2 to c3 do
        outline := outline + Format('%4d', [f[i-1, j-1]]);
      outline := outline + Format('%7d', [rowtot[i-1]]);
      AReport.Add(outline);
    end;

    outline := 'Total         ';
    for j := c2 to c3 do
      outline := outline + Format('%4d', [S[j-1]]);
    outline := outline + Format('%7d', [T]);
    AReport.Add(outline);
    AReport.Add( '');

    if (c3 = C1) then
      done := true
    else begin
      c2 := c3 + 1;
      c3 := c2 + 15;
      if (c3 > C1) then c3 := C1;
    end;
  end; // end while not done
end; // end sub frequencies


procedure TRaschForm.GetLogs(const L, L1,L2, g, g2, f2: DblDyneVec;
  const rowtot: IntDyneVec; k : integer; const s5, S: IntDyneVec;
  T, r, C1 : integer; var v1, v2: Double; AReport: TStrings);
var
  tx, rowtx, rx, t2, t3, e : double;
  i, j : integer;
begin
  t2 := 0.0;
  tx := T;
  rx := r;
  for i := 0 to r-1 do
  begin
    rowtx := rowtot[i];
    L[i] := ln(rowtx / (tx - rowtx));
    t2 := t2 + L[i];
  end;
  t2 := t2 / rx;
  for i := 0 to r-1 do
  begin
    L1[i] := L[i] - t2;
    L2[i] := L1[i] * L1[i];
    v1 := v1 + L2[i];
  end;
  v1 := v1 / rx;

  AReport.Add('Item  Log Odds  Deviation  Squared Deviation');
  AReport.Add('----  --------  ---------  -----------------');
             //xxxx   xxxxxx    xxxxxx         xxxxxx
  for i := 0 to r-1 do
    AReport.Add('%4d   %6.2f    %6.2f         %6.2f', [i+1, L[i], L1[i], L2[i]]);

  t3 := 0.0;
  v2 := 0.0;
  for j := 0 to C1-1 do
  begin
    e := s5[j];
    g[j] := ln(e / (k - e));
    g2[j] := S[j] * g[j];
    t3 := t3 + g2[j];
    f2[j] := S[j] * (g[j] * g[j]);
    v2 := v2 + f2[j];
  end;
  t3 := t3 / tx;
  v2 := v2 / (tx - (t3 * t3));

  AReport.Add('');
  AReport.Add('Score  Frequency  Log Odds  Freq.x Log  Freq.x Log Odds Squared');
  AReport.Add('-----  ---------  --------  ----------  -----------------------');
             // xxx      xxx      xxxxxx     xxxxx             xxxxxx
  for j := 0 to C1-1 do
    AReport.Add(' %3d      %3d      %6.2f     %6.2f             %6.2f', [s5[j], S[j], g[j], g2[j], f2[j]]);
  AReport.Add('');
end;


procedure TRaschForm.GetScores(NoSelected: integer; const Selected: IntDyneVec;
  NoCases: integer; const f: IntDyneMat; const Mean, XSqr, SumXY: DblDyneVec;
  const S, X: IntDyneVec; out sumx, sumx2: double; var N: integer;
  AReport: TStrings);
var
  i, j, k1, T, item: integer;
  outline: string;
begin
   AReport.Add('RASCH ONE-PARAMETER LOGISTIC TEST SCALING (ITEM RESPONSE THEORY)');
   AReport.Add('Written by William G. Miller');
   AReport.Add('');

   k1 := NoSelected;
   for i := 1 to k1 do
   begin
     for j := 1 to k1 + 2 do
       f[i-1,j-1] := 0;
     Mean[i-1] := 0.0;
     XSqr[i-1] := 0.0;
     SumXY[i-1] := 0.0;
   end;
   for j := 1 to k1 + 2 do
     S[j-1] := 0;
   N := 0;
   SumX := 0.0;
   SumX2 := 0.0;

   // Read each case and scores for each item.  Eliminate rows (subjects)
   // that have a total score of zero or all items correct
   for i := 1 to NoCases do
   begin
     if not GoodRecord(OS3MainFrm.DataGrid, i, Selected) then
       continue;

     T := 0;
     for j := 1 to k1 do
     begin
       item := Selected[j-1];
       X[j-1] := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[item, i])));
       T := T + X[j-1];
     end;

     if (T < k1) and (T > 0) then
     begin
       outline := format('Case %3d:  Total Score: %3d  Item scores', [i, T]);
       SumX := SumX + T;
       SumX2 := SumX2 + T * T;
       for j := 0 to k1-1 do
       begin
         Mean[j] := Mean[j] + X[j];
         XSqr[j] := XSqr[j] + X[j] * X[j];
         SumXY[j] := SumXY[j] + X[j] * T;
         outline := outline + Format('%2d', [X[j]]);
         if (X[j] = 0) then
           f[j,T-1] := f[j,T-1] + 1;
       end;
       AReport.Add(outline);

       S[T-1] := S[T-1] + 1;
       N := N + 1;
     end else
       AReport.Add('Case %3d eliminated.  Total score was %3d', [i, T]);
   end;
   AReport.Add('');
end;


procedure TRaschForm.InBtnClick(Sender: TObject);
var
  i: integer;
begin
  i := 0;
  while i < VarList.Items.Count do
  begin
    if VarList.Selected[i] then
    begin
      ItemList.Items.Add(VarList.Items[i]);
      VarList.Items.Delete(i);
      i := 0;
    end else
      i := i + 1;
  end;
  UpdateBtnStates;
end;


procedure TRaschForm.Maxability(const expdcnt, d2, e2: DblDyneVec;
  const p1: DblDyneMat; const p2, P: DblDyneVec; C1, r: integer;
  const D: DblDyneMat; const s5: IntDyneVec; noloops: integer; AReport: TStrings);
var
  i, j: integer;
  d9: double;
begin
  d9 := 0.0;
  AReport.Add('Maximum Likelihood Iteration Number %d', [noloops]);
  for j := 0 to C1-1 do
  begin
    expdcnt[j] := 0.0;
    d2[j] := 0.0;
  end;
  for i := 0 to r-1 do
    for j := 0 to C1-1 do
      p1[i,j] := exp(p2[j] - P[i]) / (1.0 + exp(p2[j] - P[i]));
  for j := 0 to C1-1 do
    for i := 0 to r-1 do
    begin
      expdcnt[j] := expdcnt[j] + p1[i,j];
      // expected number in score group
      D[i,j] := exp(p2[j] - P[i]) / (sqrt(1.0 + exp(p2[j] - P[i])));
      d2[j] := d2[j] + D[i,j]; // rate of change value
    end;
  for j := 0 to C1-1 do
  begin
    e2[j] := expdcnt[j] - s5[j];  // discrepency
    e2[j] := e2[j] / d2[j];
    if (abs(e2[j]) > d9) then d9 := abs(e2[j]);
    p2[j] := p2[j] - e2[j];
  end;
{    Debug check in old sample program
	'     writeln
	'     writeln('Actual and Estimated Scores')
	'     writeln
	'     writeln('Score  Estimated  Adjustment')
	'     for j := 1 to c1 do
	'         writeln(s5(j):3,'   ',expdcnt(j):6:2,'       ',e2(j):6:2)
	'     writeln
}
end; // end of maxability


function TRaschForm.MaxItem(const R1, d1: DblDyneVec; const p1, D: DblDyneMat;
  const e1, p2, P: DblDyneVec; const S, rowtot: IntDyneVec;
  T, r, C1: integer): double;
var
  i, j: integer;
  d9: double;
begin
  d9 := 0.0;
  for i := 0 to r-1 do
  begin
    R1[i] := 0.0;
    d1[i] := 0.0;
  end;
  for i := 0 to r-1 do
    for j := 0 to C1-1 do
      p1[i,j] := exp(p2[j] - P[i]) / (1.0 + exp(p2[j] - P[i]));
  for i := 0 to r-1 do
  begin
    for j := 0 to C1-1 do R1[i] := R1[i] + S[j] * p1[i,j];
    e1[i] := R1[i] - (T - rowtot[i]);
  end;

  // e1(i) contains the difference between actual and expected passes
  // now calculate derivatives and adjustments
  for i := 0 to r-1 do
  begin
    for j := 0 to C1-1 do
    begin
      D[i,j] := exp(p2[j] - P[i]) / (sqrt(1.0 + exp(p2[j] - P[i])));
      d1[i] := d1[i] + (S[j] * D[i,j]);
    end;
    e1[i] := e1[i] / d1[i];

    // adjustment for item difficulty estimates
    if (abs(e1[i]) > d9) then d9 := abs(e1[i]);
    P[i] := P[i] + e1[i];
  end;

{    debug check from old sample program
	'     writeln
	'     writeln('actual and estimated items right')
	'     writeln
	'     writeln('item  actual  estimated  adjustment')
	'     for i := 1 to r do
	'     begin
	'          writeln(i:3,'  ',(t-rowtot(i)):3,'           ',e1(i):6:2)
	'     end
	'     writeln
}

  Result := d9;
end;

procedure TRaschForm.MAXOUT(r, C1: integer; const i5, s5: IntDyneVec;
  const P, p2: DblDyneVec; AReport: TStrings);
var
  i, j: integer;
begin
  AReport.Add('');
  AReport.Add('Maximum Likelihood Estimates');
  AReport.Add('');

  AReport.Add('Item   Log Difficulty');
  AReport.Add('----   --------------');
             //xxx        xxxxxx
  for i := 0 to r-1 do
    AReport.Add('%3d        %6.2f', [i5[i], P[i]]);

  AReport.Add('');
  AReport.Add('Score   Log Ability');
  AReport.Add('-----   -----------');
             // xxx       xxxxxx
  for j := 0 to C1-1 do
    AReport.Add(' %3d       %6.2f', [s5[j], p2[j]]);
end;


procedure TRaschForm.OutBtnClick(Sender: TObject);
var
  i: integer;
begin
  i := 0;
  while i < ItemList.Items.Count do
  begin
    if ItemList.Selected[i] then
    begin
      VarList.Items.Add(ItemList.Items[i]);
      ItemList.Items.Delete(i);
      i := 0;
    end else
      i := i + 1;
  end;
  UpdateBtnStates;
end;


procedure TRaschForm.Plot(const xyArray: DblDyneMat; Arraysize: integer);
var
  ser: TChartSeries;
  i: integer;
begin
  ser := FItemFuncsChartFrame.PlotXY(ptLines, nil, nil, nil, nil, '', DATA_COLORS[0]);
  for i := 0 to ArraySize-1 do
    ser.AddXY(xyArray[i, 0], xyArray[i, 1]);
  ser.Active := false;
end;


procedure TRaschForm.PlotInfo(r: integer; const Info, A: DblDyneMat;
  const Slope, P: DblDyneVec);
const
  MIN = -3.5;
  MAX = +3.5;
var
  cg, hincrement, Ymax, elg, term1, term2, jx: double;
  i, j, jj, size: integer;
  TestInfo: DblDyneMat = nil;
begin
  // Prepare charts
  FChartCombobox.Items.Clear;
  FItemFuncsChartFrame.Clear;
  FItemFuncsChartFrame.SetXTitle('log ability');
  FItemFuncsChartFrame.SetYTitle('Info');

  size := 0;
  hIncrement := (MAX - MIN) / 50;

  SetLength(TestInfo, 52,2);
  cg := 0.2;
  Ymax := 0;
  for i := 1 to r do // item loop
  begin
    TestInfo[i-1, 0] := 0.0;
    TestInfo[i-1, 1] := 0.0;
    jj := 1;
    jx := MIN;
    while (jx <= MAX + hIncrement) do
    begin
      if (slope[i-1] > 30) then slope[i-1] := 30;
      elg := 1.7 * slope[i-1] * (P[i-1] - jx);
      elg := exp(elg);
      term1 := 2.89 * sqr(slope[i-1] * (1.0 - cg));
      term2 := (cg + elg) * sqr(1.0 + 1.0 / elg);
      info[i-1, jj-1] := term1 / term2;
      if (info[i-1, jj-1] > Ymax) then
        Ymax := info[i-1, jj-1];
      jj := jj + 1;
      jx := jx + hIncrement;
    end;
    size := jj-1;
  end;

  // Item loop
  for i := 0 to r-1 do
  begin
    for j := 0 to size-1 do
    begin
      A[j, 0] := MIN + (hIncrement * (j+1));
      A[j, 1] := info[i, j];
      TestInfo[j, 1] := TestInfo[j, 1] + info[i, j];
    end;
    if ItemFuncsChk.Checked then
    begin
      FChartCombobox.Items.Add(Format('Item %d', [i+1]));
      Plot(A, size);
    end;
  end;

  for j := 0 to size-1 do
    TestInfo[j, 0] := MIN + hIncrement * (j+1);

  if TestInfoChk.Checked then
    PlotTest(TestInfo, size, 'Item Information Function for Test');

  FChartCombobox.ItemIndex := 0;
  SelectItemFuncsPlot(FChartCombobox);
  FItemFuncsChartFrame.Chart.Legend.Visible := false;
end;


procedure TRaschForm.PlotItems(r: integer; const i5: IntDyneVec;
  const P: DblDyneVec);
var
  i: integer;
  ser: TChartSeries;
begin
  FItemDiffsChartFrame.Clear;
  FItemDiffsChartFrame.SetTitle('LOG DIFFICULTIES FOR ITEMS');
  FItemDiffsChartFrame.SetXTitle('Item');
  FItemDiffsChartFrame.SetYTitle('Log Difficulty');
  FItemDiffsChartFrame.HorLine(0, clBlack, psSolid, '');

  ser := FItemDiffsChartFrame.PlotXY(ptBars, nil, nil, nil, nil, '', DATA_COLORS[0]);
  for i := 0 to r-1 do
    ser.AddXY(i5[i], P[i], '', DATA_COLORS[i mod Length(DATA_COLORS)]);

  FItemDiffsChartFrame.Chart.BottomAxis.Marks.Source := ser.Source;
  FItemDiffsChartFrame.Chart.BottomAxis.Marks.Style := smsXValue;
  FItemDiffsChartFrame.Chart.Legend.Visible := false;
end;


procedure TRaschForm.PlotScrs(C1: integer; const s5: IntDyneVec; const
  p2: DblDyneVec);
var
  i: integer;
  ser: TChartSeries;
begin
  FScoresChartFrame.Clear;
  FScoresChartFrame.SetTitle('LOG ABILITIES FOR SCORE GROUPS');
  FScoresChartFrame.SetXTitle('Score');
  FScoresChartFrame.SetYTitle('Log Ability');
  FScoresChartFrame.HorLine(0, clBlack, psSolid, '');
  ser := FScoresChartFrame.PlotXY(ptBars, nil, nil, nil, nil, '', DATA_COLORS[0]);
  for i := 0 to C1-1 do
    ser.AddXY(s5[i], p2[i], '', DATA_COLORS[i mod Length(DATA_COLORS)]);
  FScoresChartFrame.Chart.BottomAxis.Marks.Source := ser.Source;
  FScoresChartFrame.Chart.BottomAxis.Marks.Style := smsXValue;
  FScoresChartFrame.Chart.Legend.Visible := false;
end;


procedure TRaschForm.PlotTest(const TestInfo: DblDyneMat;
  ArraySize: integer; const Title: string{; Vdivisions, Hdivisions: integer});
var
  i: integer;
  ser: TChartSeries;
begin
  FTestInfoChartFrame.Clear;
  FTestInfoChartFrame.SetTitle(Title);
  FTestInfoChartFrame.SetXTitle('log ability');
  FTestInfoChartFrame.SetYTitle('Info');

  ser := FTestInfoChartFrame.PlotXY(ptLines, nil, nil, nil, nil, '', DATA_COLORS[0]);
  for i := 0 to ArraySize-1 do
    ser.AddXY(TestInfo[i, 0], TestInfo[i, 1]);

  FTestInfoChartFrame.Chart.Legend.Visible := false;
end;


procedure TRaschForm.Prox(const P, p2: DblDyneVec; k, r, C1 : integer;
  const L1: DblDyneVec; yexpand, xexpand: double; const g: DblDyneVec;
  T: integer; const rowtot, i5, s5: IntDyneVec; AReport: TStrings);
var
  tx, rowtx, errorterm, stdError: double;
  i, j: integer;
begin
  for i := 0 to r-1 do P[i] := L1[i] * yexpand;
  for j := 0 to C1-1 do p2[j] := g[j] * xexpand;

  if not ProxChk.Checked then
    exit;

  AReport.Add('');
  AReport.Add('Prox values and Standard Errors' );
  AReport.Add('');
  AReport.Add('Item   Scale Value   Standard Error');
  AReport.Add('----   -----------   --------------');
             //xxx      xxxxxxx        xxxxxxx

  tx := T;
  for i := 0 to r-1 do
  begin
    rowtx := rowtot[i];
    errorterm := tx / ((tx - rowtx) * rowtx);
    stdError := yexpand * sqrt(errorterm);
    if ProxChk.checked then
      AReport.Add('%3d      %7.3f        %7.3f', [i5[i], P[i], stdError]);
  end;
  AReport.Add('Y expansion factor: %8.4f', [yexpand]);
  AReport.Add('');
  AReport.Add('Score   Scale Value   Standard Error');
  AReport.Add('-----   -----------   --------------');
            //  xxx      xxxxxxx         xxxxxxx

  for j := 0 to C1-1 do
   begin
     stdError := xexpand * sqrt(k / (s5[j] * (k - s5[j])));
     AReport.Add(' %3d      %7.3f         %7.3f', [s5[j], p2[j], stdError]);
   end;
  AReport.Add('X expansion factor: %8.4f', [xexpand]);
  AReport.Add('');
end;


function TRaschForm.Reduce(k: integer; out r, T, C1: Integer;
  const i5, RowTot, s5: IntDyneVec; const f: IntDyneMat; const S: IntDyneVec;
  AReport: TStrings): integer;
var
  done: boolean;
  check, i, j, column, row: integer;
begin
  // NOW REDUCE THE MATRIX BY ELIMINATING 0 OR 1 ROWS AND COLUMNS
  AReport.Add('');

  //Store item numbers in i5 array and initialize row totals
  for i := 0 to k-1 do
  begin
    i5[i] := i+1;
    rowtot[i] := 0;
  end;

  //Store group numbers in s5 array
  r := k;
  T := 0;
  C1 := k - 1; // No. of score groups (all correct group eliminated)
  for j := 0 to C1-1 do
  begin
    s5[j] := j+1;
    T := T + S[j];
  end;

  //Get row totals of the failures matrix (item totals)
  for i := 0 to r-1 do
    for j := 0 to C1-1 do rowtot[i] := rowtot[i] + f[i,j];

  // now check for item elimination
  done := false;
  while not done do
  begin
    for i := 1 to r do
    begin
      if (rowtot[i-1] = 0) or (rowtot[i-1] = T) then
      begin
        AReport.Add('Row %3d for item %3d eliminated.', [i, i5[i-1]]);
        if (i < r) then
        begin
          for j := i to r-1 do //move rows up to replace row i
          begin
            for column := 1 to C1 do
              f[j-1, column-1] := f[j, column-1];
              rowtot[j-1] := rowtot[j];
              i5[j-1] := i5[j];
          end;
        end;
        r := r - 1;
      end; // end if
    end; // end for i

    check := 1;
    for i := 0 to r-1 do
      if ((rowtot[i] = 0) or (rowtot[i] = T)) then check := 0;
    if (check = 1) then
      done := true;
  end;

  // check for group elimination
  done := false;
  j := 1;
  while (not done) do
  begin
    if (S[j-1] = 0) then
    begin
      AReport.Add('Column %3d score group %3d eliminated - total group count = %3d', [j, s5[j-1], S[j-1]]);
      if (j < C1) then
      begin
        for i := j to C1 - 1 do
        begin
          for row := 1 to r do
            f[row-1, i-1] := f[row-1, i];
          S[i-1] := S[i];
          s5[i-1] := s5[i];
        end;
        C1 := C1 - 1;
      end
      else
        C1 := C1 - 1;
    end;

    if C1 = 0 then
    begin
      MessageDlg('Too many cases or variables eliminated', mtError, [mbOK], 0);
      FReportFrame.DisplayReport(AReport);
      AReport.clear;
      Result := 1;
      exit;
    end;

    if (S[j-1] > 0) then
      j := j + 1;

    if (j >= C1) then
    begin
      while (S[C1-1] <= 0) do
      begin
        C1 := C1 - 1;
        if C1 = 0 then
        begin
          MessageDlg('Too many cases or variables eliminated', mtError, [mbOK], 0);
          FReportFrame.DisplayReport(AReport);
          AReport.Clear;
          Result := 1;
          exit;
        end;
      end;
      done := true;
    end;
  end;

  AReport.Add('Total number of score groups: %4d', [C1]);
  AReport.Add('');
  Result := 0;
end; // end of reduce


procedure TRaschForm.Reset;
var
  i: integer;
begin
  inherited;

  if FItemDiffsChartFrame <> nil then FItemDiffsChartFrame.Clear;
  if FScoresChartFrame <> nil then FScoresChartFrame.Clear;
  if FItemFuncsChartFrame <> nil then FItemFuncsChartFrame.Clear;
  if FTestInfoChartFrame <> nil then FTestInfoChartFrame.Clear;

  VarList.Clear;
  for i := 1 to NoVariables do
    VarList.Items.Add(OS3MainFrm.DataGrid.Cells[i,0]);
  ItemList.Clear;

  ProxChk.Checked := false;
  PlotItemDiffChk.Checked := false;
  PlotScrsChk.Checked := false;
  ItemFuncsChk.Checked := false;
  TestInfoChk.Checked := false;

  UpdateBtnStates;
end;


procedure TRaschForm.SelectItemFuncsPlot(Sender: TObject);
var
  i, index: Integer;
begin
  index := FChartCombobox.ItemIndex;
  for i := 0 to FItemFuncsChartFrame.Chart.SeriesCount-1 do
    FItemFuncsChartFrame.Chart.Series[i].Active := (i = index);
  FItemFuncsChartFrame.SetTitle(Format('Item Information Function for Item %d', [index+1]));
end;


procedure TRaschForm.Slopes(const rptbis, rbis, slope: DblDyneVec;
  N: integer; sumx, sumx2: double; const sumxy: DblDyneVec; r: integer;
  const xsqr, mean: DblDyneVec);
var
  propi, term1, term2, z, Y: double;
  j: integer;
begin
  z := 0.0;
  term1 := N * sumx2 - sumx * sumx;
  for j := 0 to r-1 do
  begin
    rptbis[j] := N * sumxy[j] - mean[j] * sumx;
    term2 := N * xsqr[j] - (mean[j] * mean[j]);
    if ((term1 > 0) and (term2 > 0)) then
      rptbis[j] := rptbis[j] / sqrt(term1 * term2)
    else
      rptbis[j] := 1.0;

    propi := mean[j] / N;
    if ((propi > 0.0) and (propi < 1.0)) then
      z := inversez(propi);
    if (propi <= 0.0) then
      z := -3.0;
    if (propi >= 1.0) then
      z := 3.0;

    Y := ordinate(z);
    if (Y > 0) then
      rbis[j] := rptbis[j] * (sqrt(propi * (1.0 - propi)) / Y)
    else
      rbis[j] := 1.0;
    if (rbis[j] <= -1.0) then
      rbis[j] := -0.99999;
    if (rbis[j] >= 1.0) then
      rbis[j] := 0.99999;
    slope[j] := rbis[j] / sqrt(1.0 - (rbis[j] * rbis[j]));
  end;
end; // end of slopes procedure


procedure TRaschForm.TestFit(r, C1: integer; const f: IntDyneMat;
  const S: IntDyneVec; const P, P2: DblDyneVec; T: integer; AReport: TStrings);
var
  ct, ch, prob: double;
  i, j: integer;
  outline: string;
begin
  AReport.Add('');
  AReport.Add('Goodness of Fit Test for Each Item');
  AReport.Add('');
  AReport.Add('Item  Chi-Squared  Degrees of  Probability');
  AReport.Add('No.   Value        Freedom     of Larger Value');
  AReport.Add('----  -----------  ----------  ---------------');
             //xxxx    xxxxxxx       xxxx         xxxxxx
  ct := 0.0;
  for i := 0 to r-1 do
  begin
    ch := 0.0;
    for j := 0 to C1-1 do
      ch := ch + (exp(p2[j] - P[i]) * f[i,j]) + (exp(P[i] - p2[j]) * (S[j] - f[i,j]));
    prob := 1.0 - chisquaredprob(ch, T - C1);
    outline := format('%4d    %7.2f       %4d         %6.4f', [i+1,ch,(T-C1),prob]);
    AReport.Add(outline);
    ct := ct + ch;
  end;
  AReport.Add('');
end;


procedure TRaschForm.UpdateBtnStates;
begin
  inherited;

  if FItemDiffsChartFrame <> nil then FItemDiffsChartFrame.UpdateBtnStates;
  if FScoresChartFrame <> nil then FScoresChartFrame.UpdateBtnStates;
  if FItemFuncsChartFrame <> nil then FItemFuncsChartFrame.UpdateBtnStates;
  if FTestInfoChartFrame <> nil then FTestInfoChartFrame.UpdateBtnStates;

  InBtn.Enabled := AnySelected(VarList);
  OutBtn.Enabled := AnySelected(ItemList);
end;


function TRaschForm.Validate(out AMsg: String; out AControl: TWinControl): Boolean;
begin
  Result := false;

  if (NoVariables < 1) then
  begin
    AMsg := 'You must have data in your data grid.';
    AControl := VarList;
    exit;
  end;

  Result := True;
end;


procedure TRaschForm.VarListSelectionChange(Sender: TObject; User: boolean);
begin
  UpdateBtnStates;
end;


end.