program simpvt c c Mark Vangel, April, 1995 c c Program `simpvt' simulates the critical values to be c used for tolerance limit factors in the limit of very c large between-batch variability (equivalently, very small c within-batch variability). These critical values should c be used instead of the Satterthwaite approximation when c the data are very unbalanced. c c This program takes as input the same input file that c is used by the `recipe' tolerance limit program. If c the input filename for `recipe' is `name.dat' then c `simpvt' will produce an output file `name.crt' which c contains the simulated critical values. This critical c value file will be read by the `recipe' program and used c instead of the Satterthwaite approximation. c parameter (mtot=1000, mlvl=100, mbch=25, mpar=10, mpts=100, $ mpivsm=1000000) implicit double precision (a-h, o-z) character*80 rec character*30 flenme character*30 flenm2 dimension x(mlvl*mpar), xx(mtot*(mpar+mbch)), xpts(mpts*mpar), $ u1(mtot*mpar), s1(mpar), v1(mpar*mpar), $ ip(mlvl), iq(mtot), wk1(mtot), wk2(mtot), $ vals(mpivsm), w(mpar) data in, iout/10, 20/ c write (*,*) " Filename (without '.dat' extension) ?" read (*,fmt='(a20)') flenme do 1 i=20, 1, -1 if (flenme(i:i) .ne. ' ') then lstnbk = i go to 2 end if 1 continue 2 continue iend = lstnbk +4 c c -- Open the data file, and the critical value file c (if a critical value file exists) c Either a .dat or a .DAT extension is allowed. c Program quits if data file cannot be opened. lcflag = 1 flenme(1:iend) = flenme(1:lstnbk) // '.dat' open (unit=in, file=flenme(1:iend), status='old', iostat=lstat) if (lstat .ne. 0) then flenme(1:iend) = flenme(1:lstnbk) // '.DAT' open(unit=in, file=flenme(1:iend), status='old', iostat=lstat) lcflag = 2 end if if (lstat .ne. 0) then write (*,*) ' SIMPVT : Error opening data file; iostat = ', $ lstat stop end if if (lcflag .eq. 1) then flenm2(1:iend) = flenme(1:lstnbk) // '.crt' else flenm2(1:iend) = flenme(1:lstnbk) // '.CRT' end if open (unit=iout, file=flenm2(1:iend), status='new',iostat=lstat) irec = 0 c c -- Read in parameters c ntot -- Total number of observations c nlvl -- Number of levels of fixed effect c nbch -- Number of batches c npar -- Number of regression parameters c npts -- Number of points on regression surface at which c the tolerance limit is to be evaluated c cont -- Probability content for tolerance limit c (=.99 for A-Basis; .90 for B-Basis) c conf -- Confidence for lower tolerance limit c (=.95 for both A- and B-Basis) c 3 continue read (in,fmt='(a80)') rec irec = irec +1 if (rec(1:1) .eq. '#') go to 3 read (unit=rec,fmt=*,iostat=lstat) ntot, nlvl, nbch, npar, $ npts, cont, conf if (lstat .ne. 0) go to 99 c c -- Make sure that the problem size doesn't exceed the c dataspace of the program c ier = 0 if (ntot .gt. mtot) then write (*,*) ' SIMPVT : Error -- ntot exceeds mtot ' ier = 1 end if if (nlvl .gt. mlvl) then write (*,*) ' SIMPVT : Error -- nlvl exceeds mlvl ' ier = 1 end if if (nbch .gt. mbch) then write (*,*) ' SIMPVT : Error -- nbch exceeds mbch ' ier = 1 end if if (npar .gt. mpar) then write (*,*) ' SIMPVT : Error -- npar exceeds mpar ' ier = 1 end if if (npts .gt. mpts) then write (*,*) ' SIMPVT : Error -- npts exceeds mpts ' ier = 1 end if if (ier .ne. 0) then write (*,*) write (*,*) ' The values in the first non-comment ' write (*,*) ' line of your input file indicate that ' write (*,*) ' performing this analysis may exceed' write (*,*) ' the array bounds of this program. To ' write (*,*) ' increase these bounds, modify the ' write (*,*) ' PARAMETER and DIMENSION statements in ' write (*,*) ' the main program source file, then ' write (*,*) ' create a new executable program ' write (*,*) stop end if c c -- Now read in the levels of the fixed effects. All c two dimensional arrays used in this program are c packed into one-dimensional arrays by columns. do 10 i=1, nlvl 4 continue read (in,fmt='(a80)') rec irec = irec +1 if (rec(1:1) .eq. '#') go to 4 read (unit=rec,fmt=*, iostat=lstat) $ (x((j-1)*nlvl+i), j=1, npar) if (lstat .ne. 0) go to 99 10 continue c c -- Next read in the data, along with the corresponding c fixed level indicators ip(i) and batch indicators iq(i). do 20 i=1, ntot 5 continue read (in,fmt='(a80)') rec irec = irec +1 if (rec(1:1) .eq. '#') go to 5 read (unit=rec,fmt=*,iostat=lstat) ip(i), iq(i) if (lstat .ne. 0) go to 99 20 continue c c -- Read in the points on the regression surface at which c the lower tolerance limits are to be calculated do 30 i=1, npts 6 continue read (in,fmt='(a)') rec irec = irec +1 if (rec(1:1) .eq. '#') go to 6 read (unit=rec,fmt=*,iostat=lstat) $ (xpts((j-1)*npts+i), j=1, npar) if (lstat .ne. 0) go to 99 30 continue c c -- Build full data matrix from unique rows do 40 i=1, ntot do 50 j=1, npar xx((j-1)*ntot+i) = x((j-1)*nlvl+ip(i)) 50 continue 40 continue c c -- SVD of design matrix xx ijob = 21 ldu = ntot ldv = npar tol = 1.d-7 call dsvdc (xx , ntot, ntot, npar, s1, wk1, u1, ldu, v1, ldv, $ wk2, ijob, info) c c -- Rank (xx) irk = 0 do 60 i=1, npar if (abs(s1(i)) .lt. tol) go to 70 irk = irk +1 60 continue 70 continue c c -- Now, simulate the ratio quantiles for each of `npts' c vectors of independent variables 7 continue write (*,*) ' Number of simulation replicates ?' read (*,*) nrep if (nrep .gt. mpivsm) then write (*,*) $ " Number of replicates should not exceed 'mpivsm'" write (*,*) $ " (Either enter smaller value, or change the " write (*,*) $ " parameter 'mpivsm' and recompile the program)" go to 7 end if write (*,*) ' Integer seed ?' read (*,*) iseed zcont = ppnd16(cont, ifault) write (*,*) write (*,*) $' SIMPVT : One-Sided Random-Effect Regression Tolerance Limits' write (*,*) ' (Version 1.0, April 1995) ' write (*,*) write (*,*) ' Simulated critical values ' write (*,*) write (*,*) ' Number of values = ',npts write (*,*) ' Number of replicates = ',nrep write (*,*) ' Seed = ',iseed write (*,*) ' Input file = ',flenme(1:iend) write (*,*) ' Output file = ',flenm2(1:iend) write (*,*) write (iout,*) ' Data: ',flenme(1:iend), $ ' Num. Reps.: ',nrep, $ ' Seed: ',iseed do 100 i=1, npts z = rnor(iseed) call dcopy (npar, xpts(i), npts, w, 1) call simrat $ (u1,s1,v1,iq,w,nbch,ntot,npar,nrep,irk,zcont,conf, $ wk1,wk2,vals,quant) write (*,*) quant write (iout,*) quant 100 continue c go to 200 99 continue c c -- Error in input file write (*,*) ' SIMPVT : Error reading input file ' write (*,*) ' record number = ',irec write (*,*) ' io-status = ',lstat write (*,*) ' bad record = ',rec 200 continue stop end