program recipe c RECIPE, a FORTRAN program for REGression Confidence Intervals c on PErcentiles (i.e., one-sided tolernace limits) for a c large class of mixed models. c c Mark Vangel, NIST, Gaithersburg, MD 20899 c Version 1.0, April 1995 c c RECIPE is a program for determining one-sided tolerance c limits from mixed models involving two variance components. c The fixed part of the model can be arbitrary, and the random c part need not be balanced in any way, except that it is c necessary that each observation have the same variance, c and that this variance equal the sum of the two components c of variance (e.g., no interactions between fixed and random c effects are permitted). The theory is documented in c c Vangel, M. G. (1995),``One-Sided $\beta$-Content c Tolerance Limits for Mixed Models'', Proceedings c of the Section on Physical and Engineering Sciences, c American Statistical Association c c Vangel, M. G. (1995), ``ANOVA Estimates for `Partially- c Balanced' Mixed Models'', draft manuscript. c c Vangel, M. G. (1995), ``Design Allowables From Regression c Models Using Data From Several Batches'', Proceedings c of the 12th Symposium on Composite Materials: Testing c and Design, American Society for Testing and Materials, c STP ???, to appear. c c The user's manual for this program is c c Vangel, M. G. (1995), ``A User's manual for RECIPE'', c draft manuscript. c c c Required subroutines: c c LINPACK -- DSVDC, Singular value decomposition. c BLAS -- Various level 1, 2, and 3 routines. BLAS routines c have been used as much as possible, for both c readibility and efficiency. c SLATEC -- DSORT, sort algorithm, c RNOR, normal pseudo-random variables c (and various low-level utility routines) c ZERO1 -- Root finder for noncentral-t inverse. ACM algorithm c which uses bisection and rational interpolation is c used. c TNC -- Lenth's algorithm for noncentral-t cdf, along with c its required subroutines. c parameter (mtot=100, mlvl=100, mbch=25, mpar=10, mpts=100) c c MLIM -- Set mlim = max(mpar*mpar, mbch*mbch, mpar+mbch) c (Some PC compilers don't allow expressions in c parameter statements; you might want to try the c statement ... c c parameter(mlim=max(mpar*mpar, mbch*mbch, mpar+mbch)) c c ... and see if your compiler accepts it) c parameter(mlim=625) c c Set the values in the parameter statement to suit the size of c your applications, or the memory constraints of your computer: c c mtot -- Maximum number of response values (Ys) c mlvl -- Maximum number of fixed effect levels (distinct c sets of Xs) c mbch -- Maximum number of batches (random effect levels) c mpar -- Maximum number of fixed parameters (regression c coefficients) c mpts -- Maximum number of points at which tolerance limit c is to be calculated c implicit double precision (a-h, o-z) logical satt character*80 rec character*30 flenme character*30 flenm2 c dimension x(1000), xx(3500), xpts(1000), $ xtx(100), xtxi(100), $ u1(1000), s1(10), v1(100), $ u2(3500), s2(35), $ v2(1225), $ xn(2500), h(10000), $ ip(100), iq(100), y(100), coef(10), $ eta0(100), eta1(100), tlm0(100), tlm1(100), $ xm(100), t(100), $ wk1(3500), wk2(625), wk3(135) c c -- If your compiler will allow it, PLEASE substitute the c following dimension statement for the above dimension c statement. If you must use the above dimension statemnt c and you change any of the parameter values refer to the c following statement and CAREFULLY adjust the above c dimension statement accordingly. I apologize for the c inelegance of this approach, but portability to c bush-league PC compilers makes it necessary. c c dimension x(mlvl*mpar), xx(mtot*(mpar+mbch)), xpts(mpts*mpar), c $ xtx(mpar*mpar), xtxi(mpar*mpar), c $ u1(mtot*mpar), s1(mpar), v1(mpar*mpar), c $ u2(mtot*(mpar+mbch)), s2(mpar+mbch), c $ v2((mpar+mbch)*(mpar+mbch)), c $ xn(mlvl*mbch), h(mlvl*mlvl), c $ ip(mtot), iq(mtot), y(mtot), coef(mpar), c $ eta0(mpts), eta1(mpts), tlm0(mpts), tlm1(mpts), c $ xm(mpts), t(mpts), c $ wk1(mtot*(mpar+mbch)), wk2(mlim), wk3(mtot+mpar+mbch), c $ simvec(mpivsm), crit(mpts) data one, zero/1.d0, 0.d0/, in1 /10/, in2 /20/ c c -- Open the data file 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) flenme((lstnbk+1):iend) = '.dat' open (unit=in1, file=flenme(1:iend), iostat=lstat, $ status='old') if (lstat .ne. 0) then flenme((lstnbk+1):iend) = '.DAT' open (unit=in1, file=flenme(1:iend), iostat=lstat, $ status='old') end if if (lstat .ne. 0) then write (*,*) 'RECIPE : Data file not found' stop end if c c -- Now check for the critical value file satt = .TRUE. flenm2 = 'n' flenm2(1:lstnbk) = flenme(1:lstnbk) flenm2((lstnbk+1):iend) = '.crt' open (unit=in2, file=flenm2(1:iend), iostat=lstat, $ status='old') if (lstat .eq. 0) satt = .FALSE. if (lstat .ne. 0) then flenm2((lstnbk+1):iend) = '.CRT' open (unit=in2, file=flenm2(1:iend), iostat=lstat, $ status='old') if (lstat .eq. 0) satt = .FALSE. end if c c -- Print out header, and indicate whether Satterthwaite c approximation is being used. write (*,*) write (*,*) $' RECIPE : One-Sided Random-Effect Regression Tolerance Limits' write (*,*) ' (Version 1.0, April 1995) ' write (*,*) if (.NOT. satt) then write (*,*) $ ' *** Simulated pivot critical values from file ', $ flenm2(1:iend),' will be used.' else write (*,*) $ ' *** Simulated pivot critical value file ', $ flenm2(1:iend),' not found.' write (*,*) ' Satterthwaite approximation will be used.' end if 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 (in1,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 (*,*) ' RECIPE : Error -- ntot exceeds mtot ' ier = 1 end if if (nlvl .gt. mlvl) then write (*,*) ' RECIPE : Error -- nlvl exceeds mlvl ' ier = 1 end if if (nbch .gt. mbch) then write (*,*) ' RECIPE : Error -- nbch exceeds mbch ' ier = 1 end if if (npar .gt. mpar) then write (*,*) ' RECIPE : Error -- npar exceeds mpar ' ier = 1 end if if (npts .gt. mpts) then write (*,*) ' RECIPE : 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 (in1,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 y(i), along with the corresponding c fixed level indicators ip(i) and batch indicators iq(i). do 20 i=1, ntot 5 continue read (in1,fmt='(a80)') rec irec = irec +1 if (rec(1:1) .eq. '#') go to 5 read (unit=rec,fmt=*,iostat=lstat) ip(i), iq(i), y(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 (in1,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 -- Calculated the tolerance limits. The first call is c to a routine which does all of the calculations c which do not depend on y. The second call c performs the remaining calculations, which do c depend on y, and returns the result in the c vector 'xm' of point estimates, and the c vector 't' of tolerance limits. call regini (nlvl, npar, ntot, nbch, npts, x, xpts, ip, iq, $ cont, conf, xx, xtx, xtxi, xn, h, u1, $ s1, v1, u2, s2, v2, tlm0, tlm1, eta0, eta1, $ satt, in2, wk1, wk2, wk3) call regdat (npar, ntot, nbch, npts, xpts, y, coef, u1, s1, $ v1, u2, tlm0, tlm1, eta0, eta1, wk1, xm, t) c write (*,*) $ ' Probability Confidence Regression Tolerance', $ ' Limit' do 40 i=1, npts write (*,1000) cont, conf, xm(i), t(i) 1000 format (2f15.2,2F15.6) 40 continue go to 100 99 continue c c -- Error in input file write (*,*) ' RECIPE : Error reading input file ' write (*,*) ' record number = ',irec write (*,*) ' io-status = ',lstat write (*,*) ' bad record = ',rec 100 continue stop end