home *** CD-ROM | disk | FTP | other *** search
/ Current Shareware 1994 January / SHAR194.ISO / engineer / ksprob21.zip / KSMISC.EXE / KSGOF2.PAS < prev    next >
Pascal/Delphi Source File  |  1993-05-19  |  5KB  |  160 lines
  1. program ksgof2; {simulates the chi-square goodness of fit test
  2.                  for a Poisson distribution with µ estimated
  3.                  from the sample}
  4.  
  5. {X is a random variable with µ. A sample of size n is taken from
  6.  the distribution of X. This program simulates the chi-square
  7.  goodness of fit test for
  8.          H0: The distribution of X is Poisson
  9.  against H1: The distribution of X is not Poisson.
  10.  
  11.  The input file must be an ascii file with the following
  12.  information in successive lines:
  13.  
  14.  1. The sample mean.
  15.  
  16.  2. The number of cells.
  17.  
  18.  3. one line for each cell, containing
  19.  
  20.     for all but the last cell, the largest value of x included
  21.     in the cell and the observed cell count.
  22.  
  23.     for the last cell, the smallest value of x included in the
  24.     cell and the observed cell count.
  25.  
  26.  When there are two numbers in one line, separate them by spaces.
  27.  
  28.  This cell structure allows adjacent values of x to be combined
  29.  into one cell. Since the sample mean cannot be computed from
  30.  cell information, it must be supplied separately.
  31.  
  32.  ksgof2.dat is a sample data set.
  33.  
  34.  The results written to screen are:
  35.  
  36.    the number of tables simulated with chi-square greater than
  37.    the chi-square for the observed data.
  38.  
  39.    the number of tables simulated.
  40.  
  41.    the upper and lower limits of a 95% confidence interval for
  42.    the pvalue of the chi-square goodness of fit test.
  43.  
  44.  Program limits:
  45.  
  46.    number of cells <= 100
  47.    sample size n <= 32767
  48.    number of simulated samples <= 2147483647
  49. }
  50.  
  51. {$B-,D+,E+,F+,L+,N+,R-,V-}
  52.  
  53. uses dos,crt,ksutils;
  54.  
  55. type
  56.   real = double;
  57.  
  58. var
  59.   origobs,simobs : array [1..100] of word;
  60.   origexp,simexp : array [1..100] of real;
  61.   maxx : array [1..100] of word;
  62.   n,ncells : word;
  63.   numlarger,ntosim : longint;
  64.   origchisq,simchisq,critz,origxbar,simxbar : real;
  65.   ifiln : string[80];
  66.   ifil : text;
  67.  
  68. procedure getdata;
  69. var
  70.   i : word;
  71. begin
  72.   clrscr;
  73.   gotoxy(1,1); write('input file: ');
  74.   readln(ifiln);
  75.   assign(ifil,ifiln); {$I-} reset(ifil); {$I+}
  76.   if ioresult <> 0 then halt;
  77.   readln(ifil,origxbar);
  78.   readln(ifil,ncells);
  79.   if ncells > 100 then begin close(ifil); halt end;
  80.   for i := 1 to ncells do readln(ifil,maxx[i],origobs[i]);
  81.   close(ifil)
  82. end; {getdata}
  83.  
  84. procedure getstarted;
  85. var
  86.   i : word;
  87. begin
  88.   clrscr;
  89.   for i := 1 to 100 do origobs[i] := 0;
  90.   for i := 1 to 100 do simobs[i] := 0;
  91.   for i := 1 to 100 do origexp[i] := 0.0;
  92.   for i := 1 to 100 do simexp[i] := 0.0;
  93.   for i := 1 to 100 do maxx[i] := 0;
  94.   getdata;
  95.   n := 0;
  96.   for i := 1 to ncells do n := n + origobs[i];
  97.   origexp[1] := n * poiscdf(maxx[1],origxbar);
  98.   for i := 2 to ncells - 1 do
  99.      origexp[i] := n * (poiscdf(maxx[i],origxbar) - poiscdf(maxx[i-1],origxbar));
  100.   origexp[ncells] := n * (1.0 - poiscdf(maxx[ncells] - 1,origxbar));
  101.   origchisq := 0.0;
  102.   for i := 1 to ncells do
  103.     origchisq := origchisq + sqr(origobs[i] - origexp[i]) / origexp[i];
  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.   i,j : word;
  115.   k : longint;
  116.   ival : integer;
  117.   ranum,cs,cumprob,halfwidth,ll,ul,pval,tmpr1,tmpr2,xbar : real;
  118. begin
  119.   for k := 1 to ntosim do
  120.   begin
  121.     for i := 1 to ncells do simobs[i] := 0;
  122.     simxbar := 0.0;
  123.     for i := 1 to n do
  124.       repeat
  125.         pval := randnum1;
  126.         ival := poisinv(pval,origxbar,tmpr1,tmpr2);
  127.         if ival <> 32767 then
  128.         begin
  129.           inc(ival);
  130.           simxbar := simxbar + ival;
  131.           j := 0;
  132.           repeat inc(j) until (ival <= maxx[j]) or (j = ncells);
  133.           inc(simobs[j])
  134.         end
  135.       until ival <> 32767;
  136.     simxbar := simxbar / n;
  137.     simexp[1] := n * poiscdf(maxx[1],simxbar);
  138.     for i := 2 to ncells - 1 do
  139.        simexp[i] := n * (poiscdf(maxx[i],simxbar) - poiscdf(maxx[i-1],simxbar));
  140.     simexp[ncells] := n * (1.0 - poiscdf(maxx[ncells] - 1,simxbar));
  141.     simchisq := 0.0;
  142.     for i := 1 to ncells do
  143.       simchisq := simchisq + sqr(simobs[i] - simexp[i]) / simexp[i];
  144.     if simchisq >= origchisq then inc(numlarger);
  145.     gotoxy(1,1); write(' ':14);
  146.     gotoxy(1,1); write(numlarger, ' ',k)
  147.   end;
  148.   halfwidth := sqrt((numlarger/k)*((k - numlarger)/k)/k)*critz;
  149.   ul := numlarger / k + halfwidth;
  150.   ll := numlarger / k - halfwidth;
  151.   clrscr; gotoxy(1,1); write(numlarger, ' ',k,' ',ll:3:10,' ',ul:3:10)
  152. end; {generatesamples}
  153.  
  154. begin {main}
  155.   clrscr;
  156.   initrand1;
  157.   getstarted;
  158.   generatesamples
  159. end.
  160.