home *** CD-ROM | disk | FTP | other *** search
/ Current Shareware 1994 January / SHAR194.ISO / engineer / ksprob21.zip / KSMISC.EXE / KSGOF3.PAS < prev    next >
Pascal/Delphi Source File  |  1993-05-19  |  4KB  |  159 lines
  1. program ksgof3; {simulation of the chi-square test for independence
  2.                  in an r x c contingency table}
  3.  
  4. {This program simulates the chi-square goodness of fit test for
  5.          H0: The row and column variables are independent
  6.  against H1: These variables are not independent.
  7.  
  8.  The input file must be an ascii file with the following
  9.  information in successive lines:
  10.  
  11.  1. The number of rows and the number of columns in the table.
  12.  
  13.  2. one line for each row, containing the observed counts
  14.     for that row.
  15.  
  16.  separate numbers in the same line with spaces.
  17.  
  18.  ksgof3.pas is a sample data set.
  19.  
  20.  The results written to screen are:
  21.  
  22.    the number of tables simulated with chi-square greater than
  23.    the chi-square for the observed data.
  24.  
  25.    the number of tables simulated.
  26.  
  27.    the upper and lower limits of a 95% confidence interval for
  28.    the pvalue of the chi-square goodness of fit test.
  29.  
  30.  Program limits:
  31.  
  32.    number of rows > 1
  33.    number of columns > 1
  34.    number of rows * number of columns <= 100
  35.    sample size n <= 32767
  36.    number of simulated samples <= 2147483647
  37. }
  38.  
  39. {$B-,D+,E+,F+,L+,N+,R-,V-}
  40.  
  41. uses dos,crt,ksutils;
  42.  
  43. type
  44.   real = double;
  45.  
  46. var
  47.   origobs,simobs : array [1..100] of word;
  48.   e : array [1..100] of real;
  49.   rsum,csum : array [1..50] of word;
  50.   n,r,c : word;
  51.   numlarger,ntosim : longint;
  52.   origchisq,simchisq,critz : real;
  53.   ifiln : string[80];
  54.   ifil : text;
  55.  
  56. procedure getdata;
  57. var
  58.   i,j : word;
  59. begin
  60.   clrscr;
  61.   gotoxy(1,1); write('input file: ');
  62.   readln(ifiln);
  63.   assign(ifil,ifiln); {$I-} reset(ifil); {$I+}
  64.   if ioresult <> 0 then halt;
  65.   readln(ifil,r,c);
  66.   if (r * c > 100) or (r < 2) or (c < 2)
  67.     then begin close(ifil); halt end;
  68.   for i := 1 to r do
  69.   begin
  70.     for j := 1 to c do read(ifil,origobs[(i - 1) * c + j]);
  71.     readln(ifil)
  72.   end;
  73.   close(ifil)
  74. end; {getdata}
  75.  
  76. procedure getstarted;
  77. var
  78.   i,j,k : word;
  79. begin
  80.   clrscr;
  81.   for i := 1 to 100 do origobs[i] := 0;
  82.   for i := 1 to 100 do simobs[i] := 0;
  83.   for i := 1 to 100 do e[i] := 0.0;
  84.   for i := 1 to 50 do rsum[i] := 0;
  85.   for i := 1 to 50 do csum[i] := 0;
  86.   getdata;
  87.   for i := 1 to r do
  88.     for j := 1 to c do
  89.     begin
  90.       k := (i - 1) * c + j;
  91.       rsum[i] := rsum[i] + origobs[k];
  92.       csum[j] := csum[j] + origobs[k]
  93.     end;
  94.   n := 0;
  95.   for i := 1 to r do n := n + rsum[i];
  96.   origchisq := 0.0;
  97.   for i := 1 to r do
  98.     for j := 1 to c do
  99.     begin
  100.       k := (i - 1) * c + j;
  101.       e[k] := rsum[i] * csum[j] / n;
  102.       origchisq := origchisq + sqr(origobs[k] - e[k]) / e[k]
  103.     end;
  104.   numlarger := 0;
  105.   critz := zinv(1.0 - 0.025,10);
  106.   clrscr;
  107.   gotoxy(1,1); write('number of tables to simulate: ');
  108.   readln(ntosim);
  109.   gotoxy(1,1); write(' ':40);
  110. end; {getstarted}
  111.  
  112. procedure generatesamples;
  113. var
  114.   simrsum,simcsum : array [1..50] of word;
  115.   i,j,cellid,ntodo : word;
  116.   k : longint;
  117.   ranum,cumprob,halfwidth,ll,ul : real;
  118. begin
  119.   for k := 1 to ntosim do
  120.   begin
  121.     for i := 1 to r * c do simobs[i] := 0;
  122.     for i := 1 to r do simrsum[i] := 0;
  123.     for j := 1 to c do simcsum[j] := 0;
  124.     for ntodo := n downto 1 do
  125.     begin
  126.       ranum := randnum1;
  127.       cumprob := 0.0;
  128.       cellid := 0;
  129.       repeat
  130.         inc(cellid);
  131.         i := (cellid - 1) div c + 1;
  132.         j := (cellid - 1) mod c + 1;
  133.         cumprob := cumprob + ((rsum[i] - simrsum[i]) / ntodo) *
  134.                              ((csum[j] - simcsum[j]) / ntodo);
  135.       until (cumprob >= ranum) or (cellid = r * c);
  136.       inc(simobs[cellid]);
  137.       inc(simrsum[i]);
  138.       inc(simcsum[j]);
  139.     end;
  140.     simchisq := 0.0;
  141.     for i := 1 to r * c do
  142.       simchisq := simchisq + sqr(simobs[i] - e[i]) / e[i];
  143.     if simchisq >= origchisq then inc(numlarger);
  144.     gotoxy(1,1); write(' ':14);
  145.     gotoxy(1,1); write(numlarger, ' ',k)
  146.   end;
  147.   halfwidth := sqrt((numlarger/k)*((k - numlarger)/k)/k)*critz;
  148.   ul := numlarger / k + halfwidth;
  149.   ll := numlarger / k - halfwidth;
  150.   clrscr; gotoxy(1,1); write(numlarger, ' ',k,' ',ll:3:10,' ',ul:3:10)
  151. end; {generatesamples}
  152.  
  153. begin {main}
  154.   clrscr;
  155.   initrand1;
  156.   getstarted;
  157.   generatesamples
  158. end.
  159.