uregm.pas

心理测量统计分析中

unit URegM;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls;

type
  TForm1 = class(TForm)
    MsgLabel: TLabel;
    Label2: TLabel;
    InFlNmEdit: TEdit;
    Label3: TLabel;
    OutFlNmEdit: TEdit;
    OKButton: TButton;
    ExitButton: TButton;
    OpenDialog1: TOpenDialog;
    procedure FormActivate(Sender: TObject);
    procedure ExitButtonClick(Sender: TObject);
    procedure OKButtonClick(Sender: TObject);
  private
    { Private 愰尵 }
  public
    { Public 愰尵 }
  end;

var
  Form1: TForm1;


const Dim1 = 100;
      Dim2 = 50;

type  TMat = array[1..Dim2,1..Dim2] of Extended;


implementation

{$R *.DFM}

(*      This program was written by Y.Okamoto, 1998,3    *)

type TMatGInv = TMat;

{$i GInv.inc }


type TMatDet  = TMat;

{$i Det.inc  }

procedure TForm1.FormActivate(Sender: TObject);
begin
    with OpenDialog1 do
      begin
          Title:='Input File';
          Execute;
          InFlNmEdit.Text:=FileName;

          Title:='Output File';
          FileName:='';
          Execute;
          OutFlNmEdit.Text:=FileName;
      end;

    MsgLabel.Caption:='Press OK-Button';
    OKButton.SetFocus;
end;

procedure TForm1.ExitButtonClick(Sender: TObject);
begin
            Close;
end;


var   y, x, ya : array[1..Dim1,1..Dim2] of Extended;
      IXGInvXXX : array[1..Dim1,1..Dim1] of Extended;
      b, XX, InvXX, h, xy, c, a, cxxc, InvCXXC,
      Qea, Qta, cba : TMat;


procedure TForm1.OKButtonClick(Sender: TObject);
Label QPZ, QP, QP1;
var   InF, OutF : TextFile;
      v, ckv, VQea, VQta, Lambda, F, Chi2 : Extended;
      n, r, p, py, i, j, i1, i2, k, rank, g, u : Longint;
      s  : string;
      cc : char;
      ck : boolean;
begin
     MsgLabel.Caption:='寁嶼拞偱偡';
     AssignFile(InF,InFlNmEdit.Text); Reset(InF);
     AssignFile(OutF,OutFlNmEdit.Text); Rewrite(OutF);

     (*   僨乕僞偺撉傒崬傒   *)

     ck:=false;
     repeat
       readln(InF, s);
       if Length(s) > 0 then
         if s[1] = '/' then ck:=true;
     until ck;

     readln(InF, py, p);
     n:=0;
     while True do
       begin
           n:=n+1;
           read(InF, cc);
           if cc = '/' then goto QP;
           for i:=1 to py do read(InF, y[n,i]);
           for i:=1 to p do read(InF, x[n,i]);
           readln(InF);
       end;
 QP : n:=n-1;    //  桳岠僨乕僞悢

     writeln(OutF,'擖椡僨乕僞僼傽僀儖 = ', InFlNmEdit.Text);
     writeln(OutF);
     for i:=1 to n do
       begin
         write(OutF, ' y = ');
         for j:=1 to py do
           write(OutF,'  ',y[i,j]:10:3);
         write(OutF,'   x = ');
         for j:=1 to p do write(OutF,' ',
                                FloatToStrF(x[i,j],ffGeneral,7,1));
         writeln(OutF);
       end;

    (*    X'X偺寁嶼   *)

     for i:=1 to p do
       for j:=i to p do
         begin
             v:=0.0;
             for k:=1 to n do v:=v+x[k,i]*x[k,j];
             XX[i,j]:=v;
             XX[j,i]:=v;
         end;

    (*   X'X偺堦斒媡峴楍偺寁嶼    *)

     GInv( XX, c, p, p, r );

     for i:=1 to p do
       for j:=i to p do
         begin
             InvXX[i,j]:=0.5*(c[i,j]+c[j,i]);
             InvXX[j,i]:=InvXX[i,j];
         end;

     for i:=1 to p do
       for j:=1 to p do
         begin
           v:=0.0;
           for k:=1 to k do
             v:=v+InvXX[i,k]*XX[k,j];
           h[i,j]:=v;
         end;

    (*   b = (GInv(X'X))X'y 偺寁嶼   *)

     for i:=1 to p do
       for j:=1 to py do
         begin
           v:=0.0;
           for i1:=1 to n do v:=v+x[i1,i]*y[i1,j];
           xy[i,j]:=v;
         end;

     for i:=1 to p do
       for j:=1 to py do
         begin
           v:=0.0;
           for i1:=1 to p do v:=v+InvXX[i,i1]*xy[i1,j];
           b[i,j]:=v;
         end;

     writeln(OutF);
     writeln(OutF, 'B =');
     for i:=1 to p do
       begin
         for j:=1 to py do
           write(OutF,'      ',
                      FloatToStrF(b[i,j],ffGeneral,9,4));
           writeln(OutF);           
       end;

     (*  I - X*GInv(XX)*X 偺寁嶼  *)

     MsgLabel.Caption:='Calculation of I-XGInvXXX started';
     UpDate;
     for i:=1 to n do
       for j:=i to n do
         begin
           v:=0.0;
           for i1:=1 to p do
             for i2:=1 to p do
               v:=v+x[i,i1]*InvXX[i1,i2]*x[j,i2];
           if i=j then IXGinvXXX[i,j]:=1.0-v
                  else
                    begin
                      IXGinvXXX[i,j]:=-v;
                      IXGinvXXX[j,i]:=-v;
                    end;
         end;
     MsgLabel.Caption:='Calculation of I-XGInvXXX ended';
     UpDate;


    (*   僐儞僩儔僗僩 CBA = 0 偺専掕   *)

     while True do
       begin
           (*   僐儞僩儔僗僩 C 偺撉傒崬傒  *)

           ck:=false;
           repeat
             readln(InF, s);
             if Length(s) > 0 then
               if s[1] = '/' then ck:=true;
           until ck;

           readln(InF, g);
           if g < 1 then goto QPZ;

           if g > r then
             begin
               writeln(OutF);
               writeln(OutF,'Rank of contrast matrix C is too large');
               writeln(OutF,'Nominal rank of C = ',g);
               writeln(OutF,'Rank of C should be smaller than or equal to ',
                             r);
               goto QP1;
             end;

           for i:=1 to g do
             begin
               for j:=1 to p do read(InF,c[i,j]);
               readln(InF);
             end;

           writeln(OutF);
           writeln(OutF,'僐儞僩儔僗僩 C =');
           for i:=1 to g do
             begin
               for j:=1 to p do write(OutF,c[i,j]:7:1);
               writeln(OutF);
             end;

        (*   Check of Estimability   *)

           ckv:=0.0;
           for i:=1 to g do
             for j:=1 to p do
               begin
                 v:=0.0;
                 for k:=1 to p do
                   v:=v+c[i,k]*h[k,j];
                 ckv:=ckv+abs(v-c[i,j]);
               end;

           if ckv > 1.0E-9 then
             begin
               writeln(OutF);
               writeln(OutF,'The Contrast Matrix Is Not Estimable');

               ck:=false;
               repeat
                 readln(InF, s);
                 if Length(s) > 0 then
                   if s[1] = '/' then ck:=true;
               until ck;

               goto QP1;
             end;

        (*   僐儞僩儔僗僩 A 偺撉傒崬傒  *)

           ck:=false;
           repeat
             readln(InF, s);
             if Length(s) > 0 then
               if s[1] = '/' then ck:=true;
           until ck;

           readln(InF, u);

           if u > py then
             begin
               writeln(OutF);
               writeln(OutF,'Rank of contrast matrix A is too large');
               writeln(OutF,'Nominal rank of A = ', u);
               writeln(OutF,'Rank of A should be smaller than or equal to ',
                             py);
               goto QP1;
             end;

           for i:=1 to py do
             begin
               for j:=1 to u do read(InF,a[i,j]);
               readln(InF);
             end;

           writeln(OutF);
           writeln(OutF,'py = ',py,'    u = ',u);  
           writeln(OutF);
           writeln(OutF,'僐儞僩儔僗僩 A =');
           for i:=1 to py do
             begin
               for j:=1 to u do write(OutF,a[i,j]:7:1);
               writeln(OutF);
             end;

           Ginv( a, cxxc, py, u, rank );

           if rank < u then
             begin
               writeln(OutF);
               writeln(OutF, 'Rank of contrast matrix A is not enough !');
               writeln(OutF, 'Actual Rank = ', rank);
               writeln(OutF, 'Required Rank = ', u);
               goto QP1;
             end;

           GInv( c, cxxc, g, p, rank );

           if rank < g then
             begin
               writeln(OutF);
               writeln(OutF, 'Rank of contrast matrix C is not enough !');
               writeln(OutF, 'Actual Rank = ', rank);
               writeln(OutF, 'Required Rank = ', g);
               goto QP1;
             end;

         (*     Qea 偺寁嶼    *)

           for i:=1 to n do
             for j:=1 to u do
               begin
                 v:=0.0;
                 for i1:=1 to py do
                   v:=v+y[i,i1]*a[i1,j];
                 ya[i,j]:=v;
               end;

           MsgLabel.Caption:='Calculation of Qea started';
           UpDate;
           for i:=1 to u do
             for j:=i to u do
               begin
                 v:=0.0;
                 for i1:=1 to n do
                   for i2:=1 to n do
                   v:=v+ya[i1,i]*IXGInvXXX[i1,i2]*ya[i2,j];
                 Qea[i,j]:=v;
                 Qea[j,i]:=v;
               end;
           MsgLabel.Caption:='Calculation of Qea ended';
           UpDate;

           VQea:=Det( Qea, u, 1.0e-15 );

           writeln(OutF);
           writeln(OutF);
           writeln(OutF,'Det(Qea) = ', FloatToStrF(VQea,ffGeneral,9,4));


           (*   cxxc = c*Inv(X'X)*c' 偺寁嶼  *)

           for i:=1 to g do
             for j:=i to g do
               begin
                 v:=0.0;
                 for i1:=1 to p do
                   for i2:=1 to p do
                     v:=v+c[i,i1]*InvXX[i1,i2]*c[j,i2];
                 cxxc[i,j]:=v;
                 cxxc[j,i]:=v;
               end;

         (*   Qta = Qha + Qea 偺寁嶼   *)

           GInv( cxxc, InvCXXC, g, g, rank );

           for i:=1 to g do
             for j:=1 to u do
               begin
                 v:=0.0;
                 for i1:=1 to p do
                   for i2:=1 to py do
                     v:=v+c[i,i1]*b[i1,i2]*a[i2,j];
                  cba[i,j]:=v;
               end;

           for i:=1 to u do
             for j:=i to u do
               begin
                 v:=0.0;
                 for i1:=1 to g do
                   for i2:=1 to g do
                     v:=v+cba[i1,i]*InvCXXC[i1,i2]*cba[i2,j];
                 Qta[i,j]:=Qea[i,j]+v;
                 Qta[j,i]:=Qta[i,j];
               end;

           VQta:=Det(Qta, u, 1.0e-15);
           writeln(OutF);
           writeln(OutF,'Det(Qta) = ', FloatToStrF(VQta,ffGeneral,9,4));

           Lambda:=VQea/VQta;
           writeln(OutF);
           writeln(OutF,'Lambda = ', FloatToStrF(Lambda,ffGeneral,9,4),
                        '    df1 = ', u,
                        '    df2 = ', g,
                        '    df3 = ', n-r);

           if u = 1 then
             begin
               F:=(n-r)*((1/Lambda)-1.0)/g;
               writeln(OutF);
               writeln(OutF,'F  = ', FloatToStrF(F,ffGeneral,9,4),
                            '    df1 = ', g, '    df2 = ', n-r);
             end;

         (*    Chi-Square approximation   *)

           Chi2:=-((n-r)-0.5*(u-g+1))*Ln(Lambda);
           writeln(OutF);
           writeln(OutF,'Chi-square = ',FloatToStrF(Chi2,ffGeneral,9,4),
                        '     df = ', u*g);

         QP1 : ;
       end;

  QPZ :
     CloseFile(InF); CloseFile(OutF);
     MsgLabel.Caption:='寁嶼傪廔椆偟傑偟偨';
     ExitButton.SetFocus;
end;

end.