C DBATCH3.FOR Revised 03JAN2002 C AUTHOR: PAUL N. SOMERVILLE C This is a FORTRAN 90 program to calculate critical values for multiple test- C ing and simultaneous confidence intervals. It is an extension of the C previous QBATCH and DBATCH Fortran programs. An interactive version of C the program is DINTER3.FOR. SAS-IML versions are QBATCH4.SAS and C QINTER4.SAS. C An executable version of DBATCH3.FOR can be found at C http://pegasus.cc.ucf.edu~somervil/home.html C It should run on any PC with WINDOWS 95 or higher. C Besides being able to calculate the critical values for a large number of C specific cases, it also has the ability to calculate constants for arbitrary C sets of contrasts or hypotheses. This includes both stepup and stepdown C testing. There are no restrictions on the number of populations, the conf- C idence levels or correlations of the estimates of the population parameters. C C Another feature of the program is its ability to calculate p-values corresp- C onding to calculated or theoretical test statistic values. C Assume a vector x of dimension k is an estimate of the parameter vector mu. C Assume x has a multivariate normal distribution with unknown mean mu and var- C iance covariance matrix vc times the scalar sigma squared. vc is assumed C known, and if sigma squared is not known, it is estimated with nu degrees of C freedom in the usual manner. Two sided confidence intervals of a contrast C c'mu may be given by C [c'x - var(c'x)**.5 *d, c'x + var(c'x)**.5 *d] C where d is a function of nu, vc, k, the confidence level,1-alpha, and the set C B of contrasts. C C The program calculates the constant d for fifteen different Multiple Compar- C ison Procedures. The estimates of the elements of the population mean mu may C be correlated. In that case the variance covariance matrix will need to be C input. C The older QBATCH programs used q as the critical value. The present C program DBATCH.FOR uses d where d=q/sqrt(2). The value d is consistent C with almost all of the literature. Use of q was consistent with the Tukey C procedure. C The program can be used to get the d-values necessary for simultaneous conf- C idence intervals for an arbitrary set of one and two sided contrasts. C The program can also be used to calculate the constants necessary for one- C step and both step-down and step-up testing. C After coordinate transformations the program makes use of Monte Carlo and C quadrature. The program makes "irep" preliminary estimates each of which C uses "mocar" random directions. These are combined to produce a single value C for d and a standard error of the estimate. (If the user sets "irep" equal C to 1, the program uses an empirical formula to obtain the standard error of C estimate). C C A modification of Brent's method (page 253 Numerical Recipes,1986) is used to C solve a nonlinear equation for q, based on the Monte Carlo sampling. The C method is iterative and stops when the value of d changes less than xacc in a C following iteration. xacc is presently set at .00001. The program specifies C an interval for Brent's method to use (brakl, braku). However, the program C uses the interval only as a beginning trial interval. C Use is made of Gauss-Legendre integration in the calculation of the value of C d. In general 15 non-empty bins (nq is the number of bins used) gives accur- C acies consistent with Monte Carlo. The program is presently set to use nq = C 50. Increasing nq results in an increase in running time, but the increase C is small compared to the generation of random directions. The program C gives d values for any nonsingular variance-covariance matrix. C The following is a list of multiple comparison procedures for which the prog- C ram can calculate critical values d. The program associates a number (swtype) C with each one. C Swtype Multiple Comparison Procedure C 1 Tukey (1953) C 2 Bofinger (1985) C 3 Hsu's MCB (1981) C 4 Dunnett (one-sided) (1955) C 5 Dunnett (two-sided) (1955) C 6 Hayter (1990) C 7 Successive ordered treatments (1-sided) Liu et al (1999) C 8 Successive ordered treatments (2-sided) Liu et al (1999) C 9 Umbrella contrasts C 10 Generalized umbrella contrasts (Bretz and Hothorn - submitted) C 11 Williams type contrasts (Bretz 1999) C 12 Marcus type contrasts (Bretz 1999) C 13 Hirotsu contrasts (Hirotsu 1997,1999) C 14 McDermott and Mudholker (1993) C 15 Isotonic contrasts (Bretz 1999) C In addition the program is capable of calculating critical values for sim- C ultaneous confidence intervals for an arbitrary set of contrasts, and for C both step-down and step-up testing of arbitrary sets of hypotheses. C Swtype C 0 Arbitrary set of contrasts (also step-down constants) C 20 Modified comparisons of successive ordered ... (Liu et al 2001) C (See note at end of program commentary) C 21 Step-up constants (one-sided case) C 22 Step-up constants (two-sided case) C RUNNING THE FORTRAN 90 PROGRAM C The program requires that two files exist before the program can be run. C QCARB.IN contains the input data for the program and QCARB.OUT will contain C the output of the program. A separate set of instructions is required for C each constant that is to be calculated. C There are two options for the calculations: C Option (i) Require the program to make 'irep' independent estimates each C using 'mocar' random directions. The estimates are then averaged C and the standard deviation of the final estimate is calculated. C Suggested values for 'irep' and 'mocar' are 10 and 10000. C If 'irep' = 1, an empirical formula is used to calculate the C standard deviation. C Option (i) is recommended. C Option (ii)Require the program to use a sufficient number of random dir- C ections so that the final estimate has an estimated standard C deviation of 'givense'. Suggested sd's are .01 or .001. C QCARB.IN must contain the following: C i) two lines common to all procedures C ii) (used only if the estimates are correlated or have unequal variances) C k lines for the lower triangular portion of the variance covariance C matrix vc. One line is used for each row. C iii) (used only for procedures 4 and 5 when the estimates are correlated C or have unequal variances, and for procedures 9 to 15 and procedure 20) C iv) (use only for swtypes 0, 20, 21 or 22) C mm lines each containing a contrast from the set B C i) The two lines common to all procedures are as follows: C swcov swtype conf givense C k seed ndenom mocar irep mm C swcov This is a switch which tells the programs whether the estimates C are uncorrelated with equal variances (swcov = 1) or whether the C variance covariance or correlation matrix will be need to be in- C cluded in the input (swcov = 0) C swtype The procedure for which a value d is required. C conf The value of 1 - alpha, or overall confidence. C givense Used for option (ii) (User specifies se for d) C For option (ii) use the value 999 (say) C k Dimension of x and mu. C seed Seed for the random number generator (an integer) C mocar For option (i) number of random directions for each estimate. C For option (ii) input any integer. C ndenom Degrees of freedom associated with estimate of sigma squared. C If sigma is known, use ndenom = -1. C irep For option (i) number of preliminary estimates. C For option (ii) input any negative integer, e.g. -999. C mm Number of contrasts. Used only for swtypes 0,20,21,22. C Otherwise input any integer, e.g. 999. C iii) if((swtype=4).or.(swtype=5)).and.(swcov=0)) C line must contain number of the control population "icontrol" C if(swtype=9) C line contains the number of the "peak" population "icontrol" C if(swtype=10) C line contains "iflo" and "ighi", the presumed range of the "peak" C population C if(11.le.swtype.le.15) C line contains n1,n2,...,nk the sample sizes for populations 1 to n. C if(swtype.eq.20) C line contains: ifix, dlower C (see note at end of program commentary) C iv) For swtypes 0,20,21,22, mm contrasts follow (the set B) C One line is used for each contrast. A line consists of the coefficients C of x1, x2, ... , xk, in that order. (For swtypes 21 or 22 see special C section on step-up claculations. C Several calculations for d may be made in the same run, with the instructions C for each run being in sequence in QCARB.IN. A line with 4 or more negative C ones, for example "-1 -1 -1 -1" is used to denote the end of the run. C Comment on swtype=20: C Using swtype=20, we solve the following problem: C Given a region bounded by the mm hyperplanes below, and under the conditions C given in the common two first lines, solve for the value of d such that the C probability content of the region is "conf". C First "mm-ifix" hyperplanes:contrast.le.d*sqrt(var(contrast)). C Last "ifix" hyperplanes: contrast.le.dlower*sqrt(var(contrast)). C Swtype=20 was used to calculate the required constants for "Combined one-side C and two-sided confidence interval procedures for successive comparisons of C ordered treatment effects" Somerville, Liu, Miwa and Hayter (2001). C Swtype=20 requires two constants on the last line of the input in QCARB.IN C They are: C ifix the number of contrasts whose RHS is fixed at the value dlower C dlower In the application for which this was designed dlower is the C d-value for swtype=7 C The interval(dlower,dupper) is the "guessed" interval for d in the subrout- C ine zbrent. Since d will always be at least as large as dlower, the program C presently sets dupper=1.15*dlower. C Comments on swtypes 21 and 22: C Swtypes 21 and 22 were designed specifically for stepup testing. (See C Dunnett and Tamhane 1991 "Stepdown multiple tests for comparing treatments C with a control in unbalanced oneway layouts" Statistics in Medicine 10, C 939-947.) Swtype 21 is used when the hypotheses are 1-sided and swtype 22 C is used when all the hypotheses are 2-sided. If there are m hypotheses to be C tested, then m constants d(1) < d(2) < ... < d(m) are needed and they must be C calculated sequentially. If we have the first j of these constants then the C program enables the j+1st to be calculated. The set B contains a contrast C corresponding to each of d(1), d(2), ... d(j). C A more detailed description of the use of the program for calculating the C constants required for step-down and step-up testing is given in an un- C published manuscript "Calculation of constants for stepwise testing". C The methodology of DBATCH is described in Somerville (1997) "Multiple C testing and simultaneous confidence intervals: calculation of constants", C Computational Statistics and Data Analysis, Vol. 25, 217-233. C For a more detailed set of program instructions and examples, see the Journal C of Statistical Software article "Multiple testing and SAS-IML programs for C computation of critical values for multiple testing and simultaneous conf- C idence intervals", Somerville and Bretz (http://www.jstatsoft.org/v06/i05/ C Additional comment: For swtypes 20,21 and 22, the "structure" of the C region is dependent on the value of d. C CALCULATION OF P-VALUES: C P-values may be calculated for a given swtype by altering two items in the C input for calculating d values. Using option (i), add 100 to the value of C the appropriate swtype. For givense, substitute the calculated or theoretical C value of the test statistic. C PROGRAMMER'S NOTE: Program QBATCH.FOR was written to get q-values. Present C program DBATCH3.FOR gets d-values where d=q/sqrt(2). While program inputs C d-values, they are converted to q-values. At the conclusion of calculations C for q-values, they are converted back to d-values (swtypes.lt.20). C The binning method is not used for swtypes 21 and 22. module comm real, dimension(:,:),allocatable::vc,hm,p,c,gg,vcc real,dimension(:),allocatable::dc,wl,pl,prop,qans,znn integer,dimension(:),allocatable::ibin,ifacee,ifaceb integer,dimension(:),allocatable::indx,irank real sb,zmd2,zmnd2,zmdzn,znd2,fconst,zm,zn,cnlogn integer k,k1,ndenom,lines,mocar,irep,mm,notempty,icontrol real zero,zone,two,zoneplus,sx,ssx,zrep,coc,sd,zmean,se real brakl,braku,conf,xacc,seed,givense,zmocar real, dimension(:,:),allocatable::contr,ggfrst real, dimension(:),allocatable::step,tte real,dimension(0:6,0:2)::bbeta real qpoff,qlower,qupper,ighi,iflo integer ifix,nw,nz,oldnw,oldnz,mmq logical swpvalue data bbeta/-6.55576,-8.1885,-8.2998,-7.5002,-7.5002,-7.2686, &-8.2998,.20617,.19807,.19615,.19964,.19964,.16873,.19615,-3.8513, &-3.9409,-3.7847,-2.9872,-2.9872,-3.4967,-3.7847/ parameter (mr=714025,ia=1366,ic=150889,rm=1./mr) end module comm use comm external diffprob character(len=8)date character(len=10)time character(len=12)cpop,cconf character(len=30)stan character(len=36)corr character(len=4)is character(len=11)qv character(len=11)seeed character(len=12)pval integer nq,swcov,swtype,mm1,mm2,iseed1 C Open(7,file='spec.out',access='sequential',form='formatted', C $recl=90,status='old') Open(8,file='qcarb.out',access='sequential',form='formatted', $recl=90,status='old') Open(9,file='qcarb.in',access='sequential',form='formatted', $recl=90,status='old') C mm is number of faces or contrasts C swcov=0 requests varcov input, swtype=0 requests contrast input. 21 read(9,*)swcov,swtype, conf,givense swpvalue=.false. zero=0. zone=1. two=2. if((swcov.eq.-1).and.(swtype.eq.-1).and.(conf.eq.-1.))goto 997 if((swcov.ne.0).and.(swcov.ne.1))then write(6,*)'INPUT ERROR: swcov read as ',swcov write(8,*)'INPUT ERROR: swcov read as ',swcov go to 997 end if if((swtype.lt.zero).or.(swtype.gt.130))then write(6,*)'INPUT ERROR: swtype read as ',swtype write(8,*)'INPUT ERROR: swtype read as ',swtype go to 997 end if if((conf.le.zero).or.(conf.ge.zone))then write(6,*)'INPUT ERROR: conf read as ',conf write(8,*)'INPUT ERROR: conf read as ',conf go to 997 end if read(9,*)k,iseed1,ndenom,mocar,irep,mm if((swtype.eq.21).and.(mm.lt.2))then write(6,*)'Error: swtype and mm incompatable' write(8,*)'Error: swtype and mm incompatable' go to 997 end if if((swtype.eq.22).and.(mm.lt.2))then write(6,*)'Error: swtype and mm incompatable' write(8,*)'Error: swtype and mm incompatable' go to 997 end if if(swtype.ge.100)then write(8,*)swcov,swtype,conf,givense write(8,*)k,iseed1,ndenom,mocar,irep,mm swtype=swtype-100 swpvalue=.true. conf=.11111 C Note that value of conf is irrelevant as it is cancelled out qq=givense*sqrt(two) irep=1 if(mocar.lt.100000)mocar=100000 end if nz=iseed1 C nw=iseed2 call date_and_time(date,time) write(6,*)'date and time ',date,' ',time write(8,*)'date and time ',date,' ',time write(6,*)swcov,swtype, conf,givense write(6,*)k,iseed1,ndenom,mocar,irep,mm if((swtype.eq.0).or.((swtype.ge.20).and.(swtype.le.22))) & allocate(contr(mm,k)) if((swtype.eq.21).or.(swtype.eq.22))then allocate(indx(mm),irank(mm)) end if allocate(vc(k,k)) C If swcov.eq.0 read input variance covariance matrix C Reads lower triangular portion one row at a time. if(swcov.eq.0)then call inputvc if(vc(1,1).lt.0)then write(6,*)'Var-cov matrix has negative diagonal' write(8,*)'Var-cov matrix has negative diagonal' go to 997 end if end if icontrol=1 C next statement needed for swtype equal 4,5 or 9 in contrast routine if((((swtype.eq.4).or.(swtype.eq.5)) &.and.(swcov.eq.0)).or.(swtype.eq.9))then read(9,*)icontrol if(icontrol.lt.0)then write(6,*)'Error: control population no. cant be lt 0' write(8,*)'Error: control population no. cant be lt 0' go to 997 end if end if C If swtype.eq.10) input lower and upper values for h if(swtype.eq.10)read(9,*)iflo,ighi if((swtype.eq.10).and.(iflo.lt.0))then write(6,*)'Error: treat no. cant be neg for swtype 10' write(8,*)'Error: treat no. cant be neg for swtype 10' go to 997 end if C If swtype is between 11 and 15, input vector of sample sizes if((swtype.ge.11).and.(swtype.le.15))then allocate(znn(k)) read(9,*)(znn(i),i=1,k) end if C Input # of 'fixed' contrasts if swtype.eq.20. C Also input dlower where dlower is d-value if ifix=mm C dlower is the d-values corresponding to swtype equal to 7. if(swtype.eq.20)then read(9,*)ifix,dlower qlower=dlower*sqrt(two) qupper=1.15*qlower C The "guessed"interval containing q is (qlower,qupper) C q will always be at least as large as qlower! if(ifix.lt.0.or.(qlower.lt.0.))then write(6,*)'Error: neither ifix or dlower can be lt 0' write(6,*)'Error: neither ifix or dlower can be lt 0' go to 997 end if if(ifix.eq.mm)then write(6,*)'Input error: ifix cant equal mm' write(8,*)'Input error: ifix cant equal mm' go to 997 end if end if C Input contrasts if swtype is 0, 20, 21 or 22. C If swtype.eq.20, last ifix contrasts must be "fixed" contrasts. if((swtype.eq.0).or.(swtype.eq.20))then kmin=min(k,4) do i=1,mm read(9,*)(contr(i,j),j=1,k) C Check for possible input error xyz=-zone do j=1,kmin xyz=max(xyz,contr(i,j)) end do if(xyz.lt.zero)then write(6,*)'WARNING!!! There may be an input error.' write(6,*)'First 3 or 4 contrast coefficients are neg.' write(8,*)'WARNING!!! There may be an input error.' write(8,*)'First 3 or 4 contrast coefficients are neg.' end if end do end if C Input contrasts for swtypes 21 or 22 C Input contrasts followed by the corresponding critical value. C For last contrast use critical value for previous contrast! if((swtype.eq.21).or.(swtype.eq.22))then allocate(step(mm),tte(mm)) read(9,*)(contr(1,j),j=1,k),step(1) do i=2,mm read(9,*)(contr(i,j),j=1,k),step(i) end do end if C Make first guesses for critical values to be calculated. if((swtype.eq.22).or.(swtype.eq.21))then step(mm)=1.1*step(mm-1) end if C if((swtype.eq.21).or.(swtype.eq.22))step=step*sqrt(two) C Start here for interactive program zmocar=mocar if((swtype.ne.4).and.(swtype.ne.5).and.(swtype.ne.9))icontrol=1 if(swtype.eq.1) mm=k*(k-1) if((swtype.eq.2).or.(swtype.eq.6)) mm=k*(k-1)/2 if((swtype.eq.3).or.(swtype.eq.4)) mm=k-1 if (swtype.eq.5) mm=2*(k-1) if(swtype.eq.7)mm=k-1 if(swtype.eq.8)mm=2*(k-1) if(swtype.eq.9)then mm1=icontrol*(icontrol-1)/2 mm2=(k+1-icontrol)*(k-icontrol)/2 mm=mm1+mm2 end if if(swtype.eq.10)mm=ighi*(ighi-1)/2+(k-iflo+1)*(k-iflo)/2 if(swtype.eq.11)mm=k-1 if(swtype.eq.12)mm=k*(k-1)/2 if(swtype.eq.13)mm=k-1 if(swtype.eq.14)mm=k-1 if(swtype.eq.15)mm=2**(k-1)-1 if((swtype.gt.0).and.(swtype.lt.20))allocate(contr(mm,k)) C The next four statements set parameter values for the program C They should be adequate for almost all practical requirements. C They of course may be changed (read comment statements) brakl = .0001 braku = 99999. xacc = .00001 nq=50 cpop=' populations' cconf=' confidence.' stan='Standard error of estimate is ' corr='d-value corresponding to population' is=' is' qv='d-value is ' seeed='. Seed is ' pval=' P-VALUE is' C set miscellaneous constants k1=k-1 zm=k1 zn=ndenom zmd2=zm/two zmnd2=(zm+zn)/two zmdzn=zm/zn znd2=zn/two if(ndenom.le.0) go to 90 cnlogn=zn*alog(zn)/two 90 sx=0. ssx=0. zrep=irep allocate(hm(k1,k),c(k1,k1)) allocate(p(k,k1)) allocate(ifacee(mm),ifaceb(mm),dc(k1)) allocate(vcc(k1,k1)) C Calculate constant for recdist distribution if(ndenom.gt.0) fconst=two*beta(zmd2,znd2) C fconst=two*beta(zmd2,znd2) C put identity matrix in varcov if swcov.ne.0 if(swcov.ne.0)call varcov C Use Helmert transformation to reduce dimension to k1=k-1 call helmert vcc=matmul(matmul(hm,vc),transpose(hm)) CC write(6,*)'Helmert ',hm deallocate(vc) C Cholesky transformation of vcc=c c(transpose) call chol deallocate(vcc) C Transformation to N(0,1) coordinates is x = p w p=matmul(transpose(hm),c) deallocate(hm,c) allocate(gg(mm,k1)) allocate(wl(nq),pl(nq),prop(nq),ibin(nq)) 200 continue 210 if((swtype.gt.0).and.(swtype.lt.20))call contrast(swtype) if((swtype.ge.11).and.(swtype.le.15))deallocate(znn) C RHS=1. call calcgg if(swtype.eq.20)then gg(mm-ifix+1:mm,:)=gg(mm-ifix+1:mm,:)/qlower*sqrt(2.) allocate(ggfrst(mm-ifix,k1)) ggfrst(1:mm-ifix,:)=gg(1:mm-ifix,:)*sqrt(2.) end if 585 continue C Generate second seed needed for random normal generator nw=iseed1 do i=1,15 nw=mod(1366*nw+150889,714025) end do C Exercise random normal generator do i=1,25 rr=rnor(nz,nw) end do if(swtype.eq.20)qpoff=sphere(qlower) if((swtype.eq.21).or.(swtype.eq.22))mmq=mm-1 if((swtype.eq.21).or.(swtype.eq.22))qpoff=sphere &(sqrt(two)*step(1)) C There are 3 CASES corresponding to zrep.lt.0, zrep=1 and zrep.gt.1. C The constants bbeta have not been developed for swtype.gt.6. C A conservative approach is used, i.e. the constants for swtype.eq.0 C are used for swtype.gt.6. if(zrep.lt.zero)then C CASE 1. USER SPECIFIES DESIRED se OF ESTIMATE mocar=10000 if(swtype.lt.20)then call gauleg(zero,zone,pl,wl,nq) call freqdist(swtype,nq,qq) end if if(swtype.eq.20)then oldnz=nz oldnw=nw call gauleg(zero,sqrt(two)/qlower,pl,wl,nq) end if if((swtype.eq.21).or.(swtype.eq.22))then oldnz=nz oldnw=nw qq=zbrent(diffprob,step(mmq),step(mmq+1),xacc,swtype,nq) end if zk=k it=swtype if(swtype.gt.6)it=0 C Since givense is for d values, we multiply by sqrt(two) to get q value qgiv=givense*sqrt(two) newmocar=1000000./qgiv**2*(exp(bbeta(it,0)+bbeta(it,1) & *qq+bbeta(it,2)/zk))**2 add=newmocar-mocar if(add.gt.0.)then tt=qq mocar=max(add,256.) if(swtype.lt.20)then call gauleg(zero,zone,pl,wl,nq) call freqdist(swtype,nq,qq) end if if(swtype.eq.20)then oldnz=nz oldnw=nw call gauleg(zero,sqrt(two)/qlower,pl,wl,nq) end if if((swtype.eq.21).or.(swtype.eq.22))then oldnz=nz oldnw=nw qq=zbrent(diffprob,step(mmq),step(mmq+1),xacc,swtype,nq) end if qq=(10000*tt+mocar*qq)/(10000+mocar) mocar=mocar+10000 end if zmean=qq se=givense if(add.lt.zero)se=exp(bbeta(it,0)+bbeta(it,1)*qq+bbeta(it,2) & /zk)*10. if(swtype.gt.6)write(8,*)'standard error is conservative' end if C END OF CASE 1 C CASE 2 ONE REPLICATION EMPIRICAL FORMULA FOR se USED. if(zrep.eq.zone)then if(swtype.lt.20)then call gauleg(zero,zone,pl,wl,nq) call freqdist(swtype,nq,qq) end if if(swtype.eq.20)then oldnz=nz oldnw=nw call gauleg(zero,sqrt(two)/qlower,pl,wl,nq) end if if((swtype.eq.21).or.(swtype.eq.22))then oldnz=nz oldnw=nw end if if(swpvalue.eqv..true.)then pvalue=zone-conf+diffprob(qq,swtype,nq) write(8,9920)pval,pvalue write(6,9920)pval,pvalue 9920 format(a12,f7.4) go to 991 end if if((swtype.eq.21).or.(swtype.eq.22))then qq=zbrent(diffprob,step(mmq),step(mmq+1),xacc,swtype,nq) end if if(swtype.le.20)qq=zbrent(diffprob,2.,5.,xacc,swtype,nq) zk=k it=swtype if(it.gt.6)it=0 zmean=qq se=exp(bbeta(it,0)+bbeta(it,1)*qq+bbeta(it,2)/zk)*1000. &/sqrt(zmocar) end if C END OF CASE 2 C CASE 3 SEVERAL REPLICATIONS USED IN OBTAINING se if(zrep.gt.zone)then sx=zero ssx=zero allocate(qans(irep)) qans=zero do njj=1,irep if(swtype.lt.20)then call gauleg(zero,zone,pl,wl,nq) call freqdist(swtype,nq,qq) end if if(swtype.eq.20)then oldnz=nz oldnw=nw call gauleg(zero,sqrt(two)/qlower,pl,wl,nq) end if if((swtype.eq.21).or.(swtype.eq.22))then oldnz=nz oldnw=nw qq=zbrent(diffprob,step(mmq),step(mmq+1),xacc,swtype,nq) end if if(swtype.le.20)qq=zbrent(diffprob,2.,5.,xacc,swtype,nq) qans(njj)=qq coc=qans(1) cm=qans(njj)-coc sx=sx+cm ssx=ssx+cm*cm end do deallocate(qans) sd=(ssx-sx*sx/zrep)/(zrep-1.) se=sqrt(sd/zrep) sd=sqrt(sd) zmean=sx/zrep+coc mocar=zrep*mocar end if C END OF CASE 3 if((swtype.eq.21).or.(swtype.eq.22)) &deallocate(indx,irank,tte,step) C Convert q values to d values for swtypes.le.20 if(swtype.le.20)then zmean=zmean/sqrt(two) se=se/sqrt(two) end if 7777 if((zrep.lt.zone).and.(swtype.eq.zero))then write(8,*)'The estimated standard error may be very conservative' write(8,*)'unless the sample means are highly correlated.' end if if((swtype.eq.3).and.(icontrol.ne.1))go to 920 write(8,9912)qv,zmean,seeed,iseed1 9912 format(a12,f8.5,a11,i7) write(8,9909)k,cpop,conf,cconf 9909 format(I3,A12,2x,F4.3,A12) if(ndenom.gt.0)then write(8,*)' df for variance estimate is ',ndenom else write(8,*)' variance is assumed known ' end if write(8,*)mocar,' random directions used. ','swtype=',swtype write(8,9910)stan,se 9910 format(a30,f9.6) go to 930 920 continue write(8,*)' ' write(8,9911)corr,icontrol,is,zmean 9911 format(a36,i4,a4,f9.6) write(8,9910)stan,se write(8,*)mocar,' random directions used. ' 930 if((icontrol.ne.k).and.(swtype.eq.3).and.(swcov.eq.0))then icontrol=icontrol + 1 mocar=zmocar go to 200 else go to 991 end if 991 continue C WE ARE NOW FINISHED AND NEED TO GO BACK for new instructions do i=1,irep end do deallocate(contr,p) if(swtype.eq.20)deallocate(ggfrst) deallocate(gg,ifacee,ifaceb,dc,wl,pl,prop,ibin) call date_and_time(date,time) write(6,*)'date and time ',date,' ',time write(8,*)'date and time ',date,' ',time write(6,*)' ' write(6,*)' ' write(8,*)'--------------------' ktype=givense 9765 format(i4,i10,i4,f6.3,f9.5,f9.5) C write(7,9765)k,ktype,ndenom,conf,zmean,se GO TO 21 C Close(7) Close(8) Close(9) 997 stop end C**************************************** Subroutine pf(n1,n2,z,poff,ier) C (This is a routine to obtain probabilities for the F-distribution) use comm integer n1,n2,ier real z,poff ier=0 if(z .gt. 0.0d0) go to 5 poff=0.0d0 ier=1 return 5 if(n1.gt.0 .and. n2.gt.0) go to 10 ier=2 poff=0.0 return 10 continue an1=n1 an2=n2 a=an1*z/(an1*z+an2) a1=1.0d0-a if(a1 .lt. 0.1d-36) a1=0.1d-36 d1=an1*0.5d0 d2=an2*0.5d0 d3=d1+d2-1.0d0 r=0.d0 s1=0.d0 s2=0.d0 del=1.0d0 xm=1.0d0 xk=1.0d0 cc=0.25d0 pi=3.141592653589793d0 n=n2 C NOTE BEGINNING OF MAJOR LOOP 15 continue C TO SEE IF DEGREES OF FREEDOM ARE ODD OR EVEN m=d2 m=2*m if(m .ne. n) go to 30 n=d2-1 C IF DEGREES OF FREEDOM ARE EVEN C N=D.F./2-1 if(n .eq. 0) go to 25 do 20 i=1,n s1=del+s1*r d2=d2-1.0d0 d3=d3-1.0d0 tem=a1/d2 r=d3*tem s2=(r+tem)*s2 20 continue 25 s1=del+s1*r del=0.d0 t=-1.0d0 d3=-1.0d0 s2=a*s2 cc=cc+0.5d0 go to 45 C IF DEGREES OF FREEDOM ARE ODD C N=(D.F.-1)/2 30 n=d2 C IF DEGREES OF FREEDOM EQUAL 1. C DO NO EXIT LOOP. if(n .eq. 0)go to 40 do 35 i=1,n s1=del+s1*r d2=d2-1.0d0 d3=d3-1.0d0 tem=a1/d2 r=d3*tem s2=(r+tem)*s2 35 continue 40 s1=xk*s1 s2=xk*s2 art=sqrt(a1) xm=xm*art t=(xm-art)/a1 d3=-0.5d0 xk=2.0d0/pi cc=cc*2.0d0 45 if(cc .gt. 0.875d0) go to 50 d2=d1 d3=d2+d3 s2=s1 s1=0.0d0 a1=a if(a1 .lt. 0.1d-36) a1=0.1d-36 n=n1 go to 15 50 if(cc .lt. 1.125d0) del=4.0d0/pi*atan(t) poff=xm*(s2-s1)-del if(0.0d0 .le. poff .and. 1.0d0 .ge. poff) return if(poff .lt. 0.0d0) poff=0.0d0 if(poff .gt. 1.0d0) poff=1.0d0 ier=3 return end C**************************************** C This routine does the binning of the reciprocal distances subroutine fillbin2(recdist,nq) use comm integer i,nq real recdist do i=nq,1,-1 if(recdist .lt. pl(i)) go to 40 j=nq+1-i ibin(j)=ibin(j)+1 go to 50 40 end do 50 return end C**************************************** C This routine calculates proportions in the bins subroutine binprop(points,nq) use comm real points integer i do i=2,nq ibin(i)=ibin(i-1)+ibin(i) end do do i=1,nq prop(i)=1.-float(ibin(i))/points end do return end C**************************************** C This routine calculates the probability of correct statements C outside of region of radius q. function fringe(nq) use comm real znm1,x,xx,cd integer j zone=1. two=2. znm1=zn-1. x=0. do i=1,nq j=nq-i+1 C cc=zm*alog(zone+zn*ppl(i)*ppl(i))/two cc=zm*alog(zone+zn*pl(i)*pl(i))/two C cd=zn*alog(zone+zone/zn/ppl(i)/ppl(i))/two cd=zn*alog(zone+zone/zn/pl(i)/pl(i))/two xx=exp(-alog(pl(i))-cc-cd) C x=x+xx*wwl(i)*prop(j) x=x+xx*wl(i)*prop(j) end do fringe=x*fconst return end C******************************** Subroutine gauleg(x1,x2,x,w,n) C From "Numerical Methods" p 125. Gives constants for Gauss-Legengre. C Original range is x1 to x2. Abscissas and weights are arrays x and w C of length N. implicit double precision(a-h,o-z) integer n real x1,x2,x(n),w(n) parameter (eps=3.d-14) m=(n+1)/2 xm=0.5d0*(x2+x1) xl=0.5d0*(x2-x1) do 12 i=1,m z=cos(3.141592654d0*(i-.25d0)/(n+.5d0)) 1 continue p1=1.0d0 p2=0.0d0 do 11 j=1,n p3=p2 p2=p1 p1=((2.0d0*j-1.0d0)*z*p2-(j-1.d0)*p3)/j 11 continue pp=n*(z*p1-p2)/(z*z-1.d0) z1=z z=z1-p1/pp if(abs(z-z1).gt.eps) go to 1 x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.0d0*xl/((1.d0-z*z)*pp*pp) w(n+1-i)=w(i) 12 continue return end C************************************** C Subroutine to generate direction cosines of a random direction in k-1 C dimensions subroutine racos use comm real sum 10 sum=0. do i=1,k1 dc(i)=rnor(nz,nw) sum=sum+dc(i)*dc(i) end do sum=sqrt(sum) do i=1,k1 dc(i)=dc(i)/sum end do return end C**************************************** FUNCTION zbrent(diffprob,x1,x2,tol,swtype,nq) C zbrent has been modified so that (x1,x2) is no longer required C to bracket the root. Assumption is made that the function C diffprob is a decreasing function. INTEGER ITMAX REAL zbrent,tol,x1,x2,diffprob,EPS integer swtype EXTERNAL diffprob PARAMETER (ITMAX=100,EPS=3.e-8) INTEGER iter REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm a=x1 b=x2 C write(8,*)'BRENT x1 and x2 are ',x1,x2 fa=diffprob(a,swtype,nq) if(fa.eq.0.)then zbrent=a go to 60 end if fb=diffprob(b,swtype,nq) if(fb.eq.0.)then zbrent=b go to 60 end if C write(8,*)'BRENT f_x1 and f_x2 are ',fa,fb icount=1 5 if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.))then C Following adjustments assume diffprob(b,swtype,nq) is a C decreasing function of b. if(fa.gt.0)then if(icount.ge.15)then write(8,*)'Impossible to get large enough q-value' write(6,*)'Impossible to get large enough q-value' go to 60 end if icount=icount+1 a=max(a,b) fa=min(fa,fb) b=a+1. fb=diffprob(b,swtype,nq) go to 5 else if(icount.ge.15)then write(8,*)'Impossible to get small enough q-value' write(6,*)'Impossible to get small enough q-value' go to 60 end if icount=icount+1 a=min(a,b) b=a/2. fb=diffprob(b,swtype,nq) go to 5 end if end if c=b fc=fb do 11 iter=1,ITMAX if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then c=a fc=fa d=b-a e=d endif if(abs(fc).lt.abs(fb)) then a=b b=c c=a fa=fb fb=fc fc=fa endif tol1=2.*EPS*abs(b)+0.5*tol xm=.5*(c-b) if(abs(xm).le.tol1 .or. fb.eq.0.)then zbrent=b return endif if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then s=fb/fa if(a.eq.c) then p=2.*xm*s q=1.-s else q=fa/fc r=fb/fc p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.)) q=(q-1.)*(r-1.)*(s-1.) endif if(p.gt.0.) q=-q p=abs(p) if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d d=p/q else d=xm e=d endif else d=xm e=d endif a=b fa=fb if(abs(d) .gt. tol1) then b=b+d else b=b+sign(tol1,xm) endif fb=diffprob(b,swtype,nq) 11 continue write(6,*)'zbrent exceeding maximum iterations' write(8,*)'zbrent exceeding maximum iterations' zbrent=b 60 return END C**************************************** function gln(a) C This function gives the logarithm of the gamma function when the C argument is a positive integer or an integer divided by 2 implicit real (a-h,o-z) i1=aint(a) c2=0. c1=1. gln=0. if((a-i1).eq.0.) then i1=i1-2 else c1=-.5 gln=alog(sqrt(acos(-1.))) endif do i=1,i1 z=c1+i x=alog(z) gln=gln+x end do return end C******************************************** function beta(a,b) C This function gives the ratio of gamma(a+b) to gamma(a)*gamma(b) beta=gln(a+b)-gln(a)-gln(b) beta= exp(beta) return end C************************************************* subroutine chol C decomposes pos def matrix s into t times t transpose use comm real,dimension(k1)::da integer i,ii,j c=0. do i=1,k1 c(i,i)=1. end do do j=1,k1 da(j)=vcc(j,j) do ii=1,j-1 da(j)=da(j)-c(j,ii)*c(j,ii)*da(ii) end do do i=j+1,k1 c(i,j)=vcc(i,j) do ii=1,j-1 c(i,j)=c(i,j)-c(j,ii)*c(i,ii)*da(ii) end do c(i,j)=c(i,j)/da(j) end do end do do j=1,k1 a=sqrt(da(j)) do i=j,k1 c(i,j)=a*c(i,j) end do end do return end C***************************************************************** subroutine helmert use comm C Calculates k-1 rows of Helmert matrix do i=1,k1 zi=i do j=k-i,k hm(i,j)=1.0d0/sqrt(zi*(zi+1.0d0)) end do do j=1,k1-i hm(i,j)=0.0d0 end do hm(i,k-i)=-zi*hm(i,k-i) end do return end C******************************************************* C SUBROUTINE TO ENUMERATE ALL OF THE PAIRS OF FACES SUBROUTINE FACELIST use comm integer i,ie,ib C DIMENSION IFACEB(KK),IFACEE(KK) C COMMON/FACE/IFACEB(190),IFACEE(190),K,KK,K1,G(190,20) I=0 IE=2 30 IB=1 40 I=I+1 IFACEB(I)=IB IFACEE(I)=IE IF(IB.NE.(IE-1)) THEN IB=IB+1 go to 40 ENDIF IF(IE.EQ.K) GO TO 90 IE=IE+1 GO TO 30 90 RETURN END C**************************************** subroutine bofdist1(dis1,dis2) C This subroutine calculate the distances from the origin to the C boundaries of a convex region with origin wholly within. use comm real dis1,dis2,a integer i,j dis1=0. dis2=0. do i=1,mm a=0. do j=1,k1 a=a+dc(j)*gg(i,j) end do if(a.lt.0.) then dis1=min(a,dis1) else dis2=max(a,dis2) endif end do dis1=-dis1 return end function fringinf(nq) use comm double precision x,xx integer i,j,nq x=0d0 do i=1,nq j=nq-i+1 cc=-gln(zmd2)-(zmd2-zone)*alog(two)-(zm+zone)*alog(pl(i)) &-zone/pl(i)/pl(i)/two xx= exp(cc) x=x+xx*wl(i)*prop(j) end do fringinf=x return end subroutine varcov C inputs or constructs var cov matrix use comm do i=1,k do j=1,k vc(i,j)= zero end do vc(i,i)=zone end do return end subroutine contrast(swtype) C Modified 19 April 1995 use comm real,dimension(:,:),allocatable::co1,co2 Comment 1 - Tukey 2 - Bof 3 - MCB Comment 4 - Dunnett 5 - Dunnett (two) 6 - Hayter (need order) Comment 7 - Order 8 - Order (two) 9 - Umbrella Comment 10- Gen. Umb 11 - Williams 12 - Marcus Comment 13- Hirotsu 14 - McDermott & Mudholkar Comment 15- Isotonic integer swtype,imm,mm1,mm2 contr=zero if((swtype.eq.1).or.(swtype.eq.2).or.(swtype.eq.6)) call facelist if(swtype.eq.1) then imm=mm/2 do i=1,imm j=i+imm contr(i,ifaceb(i))=zone contr(i,ifacee(i))=-zone contr(j,ifaceb(i))=-zone contr(j,ifacee(i))=zone end do end if if((swtype.eq.2).or.(swtype.eq.6)) then do i=1,mm contr(i,ifaceb(i))=zone contr(i,ifacee(i))=-zone end do end if if((swtype.eq.3).or.(swtype.eq.4)) then do i=1,mm contr(i,icontrol)=-zone end do do i=1,icontrol-1 contr(i,i)=zone end do do i=icontrol+1,mm contr(i,i+1)=zone end do if(icontrol.le.mm) contr(icontrol,icontrol+1)=zone end if if(swtype.eq.5) then do i=1,k1 j=i+k1 contr(i,icontrol)=-zone contr(j,icontrol)=zone end do do i=1,icontrol-1 j=i+k1 contr(i,i)=zone contr(j,i)=-zone end do do i=icontrol+1,k1 j=i+k1 contr(i,i+1)=zone contr(j,i+1)=-zone end do if(icontrol.le.k1)then contr(icontrol,icontrol+1)=zone contr(icontrol+k1,icontrol+1)=-zone end if end if if((swtype.eq.7).or.(swtype.eq.8))then do i=1,k-1 contr(i,i)=zone contr(i,i+1)=-zone end do end if if(swtype.eq.8)then do i=k,mm contr(i,i-k+1)=-zone contr(i,i-k+2)=zone end do end if if(swtype.eq.9)then mm1=icontrol*(icontrol-1)/2 mm2=(k+1-icontrol)*(k-icontrol)/2 contr=zero if(mm1.eq.0)go to 110 ifaceb=zero ifacee=zero i=0 ie=2 30 ib=1 40 i=i+1 ifaceb(i)=ib ifacee(i)=ie if(ib.ne.(ie-1))then ib=ib+1 go to 40 end if if(ie.eq.icontrol)go to 90 ie=ie+1 go to 30 90 do jj=1,mm1 contr(jj,ifaceb(jj))=-zone contr(jj,ifacee(jj))=zone end do if(mm2.eq.0)go to 910 110 ifaceb=zero ifacee=zero i=0 ie=2 300 ib=1 400 i=i+1 ifaceb(i)=ib ifacee(i)=ie if(ib.ne.(ie-1))then ib=ib+1 go to 400 end if if(ie.eq.(k-icontrol+1))go to 900 ie=ie+1 go to 300 900 do jj=1,mm2 contr(jj+mm1,ifaceb(jj)+icontrol-1)=zone contr(jj+mm1,ifacee(jj)+icontrol-1)=-zone end do end if if(swtype.eq.10)then in=0 do i=1,ighi-1 do j=i+1,ighi in=in+1 contr(in,i)=zone contr(in,j)=-zone end do end do do j=iflo,k-1 do i=j+1,k in=in+1 contr(in,j)=-zone contr(in,i)=zone end do end do end if if(swtype.eq.11)then do i=1,mm contr(i,1)=-zone j=k-i+1 zsum=sum(znn(j:k)) contr(i,j:k)=znn(j:k)/zsum end do end if if(swtype.eq.12)then allocate(co1(k1,k),co2(k1,k)) co1=zero co2=zero ij=0 do i=1,k1 j=i+1 zsum1=sum(znn(j:k)) zsum2=sum(znn(1:i)) co1(i,j:k)=znn(j:k)/zsum1 co2(i,1:i)=znn(1:i)/zsum2 end do do i=1,k1 do j=1,i ij=ij+1 contr(ij,:)=co1(i,:)-co2(j,:) end do end do deallocate(co1,co2) end if if(swtype.eq.13)then do i=1,mm zsum1=sum(znn(1:i)) zsum2=sum(znn(i+1:k)) contr(i,1:i)=-znn(1:i)/zsum1 contr(i,i+1:k)=znn(i+1:k)/zsum2 end do end if if(swtype.eq.14)then do i=1,mm zsum=sum(znn(1:i)) contr(i,1:i)=-znn(1:i)/zsum contr(i,i+1)=zone end do end if if(swtype.eq.15)call iso(k,mm,znn,contr) 910 continue return end C gammp, gammq,gser,gcf and gammln are from "Numerical Methods". FUNCTION GAMMP(A,X) IF(X.LT.0..OR.A.LE.0.)write(6,*)'x.lt.0..or.a.le.0 in GAMMP' IF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMP=GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMP=1.-GAMMCF ENDIF RETURN END FUNCTION GAMMQ(A,X) IF(X.LT.0..OR.A.LE.0.)write(6,*)'x.lt.0..or.a.le.0 in GAMMQ' IF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMQ=1.-GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMQ=GAMMCF ENDIF RETURN END SUBROUTINE GSER(GAMSER,A,X,GLN) PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) IF(X.LE.0.)THEN IF(X.LT.0.)write(6,*)'x.lt.0 in GSER' GAMSER=0. RETURN ENDIF AP=A SUM=1./A DEL=SUM DO 11 N=1,ITMAX AP=AP+1. DEL=DEL*X/AP SUM=SUM+DEL IF(ABS(DEL).LT.ABS(SUM)*EPS)GO TO 1 11 CONTINUE write(6,*)'A too large, ITMAX too small' 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN) RETURN END SUBROUTINE GCF(GAMMCF,A,X,GLN) PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) GOLD=0. A0=1. A1=X B0=0. B1=1. FAC=1. DO 11 N=1,ITMAX AN=FLOAT(N) ANA=AN-A A0=(A1+A0*ANA)*FAC B0=(B1+B0*ANA)*FAC ANF=AN*FAC A1=X*A0+ANF*A1 B1=X*B0+ANF*B1 IF(A1.NE.0.)THEN FAC=1./A1 G=B1*FAC IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1 GOLD=G ENDIF 11 CONTINUE write(6,*)'A too large, ITMAX too small in GCF' 1 GAMMCF=EXP(-X+A*ALOG(X)-GLN)*G RETURN END FUNCTION GAMMLN(XX) REAL COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP/76.18009173,-86.50532033,24.01409822, * -1.231739516,.120858003E-2,-.536382E-5,2.50662827465/ DATA HALF,ONE,FPF/0.5,1.0,5.5/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+LOG(STP*SER) RETURN END **************************************************************************** function sphere(ff) C Calculates probability content of sphere of radius ff/sqrt(2). use comm xxx=zm/two if(ndenom.gt.0)then qqdk1=ff*ff/zm/two call pf(k1,ndenom,qqdk1,poff,ier) else qqdk1=ff*ff/two/two if(qqdk1.ge.xxx)then poff=gammp(xxx,qqdk1) else poff=zone-gammq(xxx,qqdk1) end if end if sphere= poff return end subroutine inputvc C Inputs lower triangular portion of covariance matrix, one C row at a time. use comm do i=1,k read(9,*)(vc(i,j),j=1,i) end do do i=1,k-1 do j=i+1,k vc(i,j)=vc(j,i) end do end do return end subroutine freqdist(swtype,nq,ff) C For a given set of contrasts, uses binning to to estimate C distribution of inverse distances from origin to boundary use comm real points integer swtype C Set weights and locations if(swtype.eq.20)call gauleg(zero,sqrt(two)/ff,pl,wl,nq) C Note that gauleg is used here only in connection with setting C up the bins and obtaining prop. C Zero bin locations ibin=0 do i=1,mocar call racos call bofdist1(dis1,dis2) call fillbin2(dis1,nq) call fillbin2(dis2,nq) end do points=mocar*2 C Calculate bin proportions call binprop(points,nq) return end c*** Monty Python Normal function rnor(z,w) integer z,w data a,b,s/1.17741,2.506628,.8857913/ integer,parameter::m1=18000,m2=65535 z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) c*** rnor uniform, -b to b, via b*.5^31*IVNI rnor=1.16724e-9*(ishft(z,16)+iand(w,m2)) if(abs(rnor) .lt. a) return z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) c*** We want y=UNI, then y1+s-2s*exp(-.5x^2), or ln(1+s-y)<-.5x^2+ln(2s), x=s*(b-x) if(alog(1.8857913-y) .lt. .5718733-.5*rnor**2) return 1 z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) v=4.656613e-10*(ishft(z,16)+iand(w,m2)) x=-alog(abs(v))/b z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) y=.5+2.328306e-10*(ishft(z,16)+iand(w,m2)) y=-alog(y) if( y+y .lt. x*x) go to 1 rnor=sign(b+x,v) return end subroutine calcgg C Obtains LHS of boundary plane equations in w coords. RHS=1. use comm gg=matmul(contr,p) do i=1,mm dtt=dot_product(gg(i,:),gg(i,:)) dtt=sqrt(dtt) gg(i,:)=gg(i,:)/dtt end do return end Function diffprob(q,swtype,nq) C Last revised 04JAN2002 C Finds the difference between "conf" and probability content of C the confidence region for the trial value of "q". C Note that for each iteration within a given replication C the beginning seeds are the same (when swtype.eq.20,21 or 22) use comm integer swtype if (swtype.eq.20)then nz=oldnz nw=oldnw gg(1:mm-ifix,:)=ggfrst(1:mm-ifix,:)/q C Zero bin locations ibin=0 do i=1,mocar call racos call bofdist1(dis1,dis2) call fillbin2(dis1,nq) call fillbin2(dis2,nq) end do points=mocar*2 C Calculate bin proportions call binprop(points,nq) qq=sqrt(two)/qlower call gauleg(zero,qq,pl,wl,nq) if(ndenom.gt.0)x=fringe(nq) if(ndenom.le.0)x=fringinf(nq) prob=qpoff+x end if if(swtype.lt.20)then qpoff=sphere(q) qq=sqrt(two)/q call gauleg(zero,qq,pl,wl,nq) if(ndenom.gt.0)x=fringe(nq) if(ndenom.le.0)x=fringinf(nq) prob=qpoff+x end if if(swtype.eq.21)then nz=oldnz nw=oldnw xxx=zm/two sum=zero step(mm)=q do i=1,mocar call racos C Calculate contrast values and order them tte=matmul(gg,dc) call indexx(mm,tte,indx) call rank(mm,indx,irank) C Find reciprocal distance to nearest hyperplane call dist212(dis1) if(dis1.eq.zero)then sum=sum+zone cycle else qsphere=zm*dis1*dis1 qsphere=zone/qsphere if(ndenom.gt.zero)then callpf(k1,ndenom,qsphere,poff,iem) sum=sum+poff else qsphere=qsphere*xxx poff=gammp(xxx,qsphere) sum=sum+poff end if end if end do points=mocar prob=sum/points end if if(swtype.eq.22)then nz=oldnz nw=oldnw xxx=zm/two sum=zero step(mm)=q do i=1,mocar call racos C Calculate contrast values and order their absolute values tte=matmul(gg,dc) C write(6,*)'tte= ',tte C 3 lines added 16DEC2001 do itte=1,mm tte(itte)=abs(tte(itte)) end do C write(6,*)'tte= ',tte call indexx(mm,tte,indx) call rank(mm,indx,irank) C Find reciprocal distance to nearest hyperplane call dist212(dis1) C write(6,*)'dis1= ',dis1 C write(6,*)k1,ndenom,qsphere,dis1,poff,zm C write(6,*)'conf= ',conf,'prob= ',prob C write(6,*)'diffprob= ',diffprob C write(6,*)'irank= ',irank C write(6,*)'indx= ',indx C write(6,*)'step= ',step C pause if(dis1.eq.zero)then sum=sum+zone cycle else qsphere=zm*dis1*dis1 qsphere=zone/qsphere if(ndenom.gt.zero)then call pf(k1,ndenom,qsphere,poff,iem) sum=sum+poff else qsphere=qsphere*xxx poff=gammp(xxx,qsphere) sum=sum+poff end if end if end do points=mocar prob=sum/points end if diffprob=conf-prob return end Subroutine iso(k,mm,znn,cm) C Generates contrasts for swtype=15 integer k,k1,anz,s,row,l,mmm,mm1,mm2 integer i,ii,j,jj,qwer real nsum,sum_n real,dimension(:,:),allocatable::mat2,mat real,dimension(:),allocatable::n1,n2,x,y,xx,yy,ssd,vecout,nn real,dimension(:),allocatable::ssd_help,y_help,n_gew,c,con real,dimension(:),allocatable::nn_help real,dimension(k)::znn real,dimension(mm,k)::cm k1=k-1 anz=2**k1 allocate (mat2(anz,k1)) allocate (mat(anz-1,k1)) cm=0. mat2=0. mat2(2,k1)=1 C Generation of all 2^k1-1 possibilities C 0: two subsequent groups are pooled C 1: subsequent groups are not pooled do i=2,k1 do ii=2**(i-1)+1,2**i mat2(ii,k1-i+1)=1 mat2(ii,k1-i+2:k1)=mat2(ii-2**(i-1),k1-i+2:k1) end do end do mat(1:anz-1,:)=mat2(2:anz,:) row=anz-1 nsum=sum(znn(2:k)) allocate (n1(k1)) allocate (n2(k1)) n1=znn(1) n2=nsum do i=2,k1 n1(i)=sum(znn(1:i)) n2(i)=sum(znn(i+1:k)) end do C Beginning of main loop: each row = one of the C 2^k1-1 possibilities is evaluated on its own do i=1,row allocate (x(k1)) allocate (y(k1)) x=1. y=mat(i,:) s=sum(mat(i,:)) allocate (ssd(s+1)) ssd=0. nsum=0 sum_n=sum(znn) allocate (nn(k1)) allocate (nn_help(k1)) nn=-y*sqrt(n1*n2)/sqrt(sum_n) nn_help=nn nn(1:size(nn))=nn_help(size(nn):1:-1) krit=0 j=1 do if(krit.ne.0)exit if (nn(j).eq.0.) then qwer=size(nn)-1 allocate (vecout(qwer)) call schnitt(nn,j,vecout,qwer) deallocate (nn) allocate (nn(qwer)) nn=vecout deallocate (vecout) else j=j+1 endif if (j.gt.size(nn)) krit=1 end do ssd(2:s+1)=nn(1:s) allocate (ssd_help(s+1)) ssd_help=ssd ssd(1:size(ssd))=ssd_help(size(ssd):1:-1) krit=0 j=1 allocate (n_gew(size(znn))) allocate (y_help(size(y))) y_help=y n_gew=znn do if(krit.ne.0)exit if (y_help(j).eq.0) then n_gew(j)=n_gew(j)+n_gew(j+1) jj=j+1 qwer=size(n_gew)-1 allocate (vecout(qwer)) call schnitt(n_gew,jj,vecout,qwer) deallocate (n_gew) allocate (n_gew(qwer)) n_gew=vecout deallocate (vecout) jj=j qwer=size(y_help)-1 allocate (vecout(qwer)) call schnitt(y_help,jj,vecout,qwer) deallocate (y_help) allocate (y_help(qwer)) y_help=vecout deallocate (vecout) else j=j+1 endif krit=1 do ii=1,size(y_help) if (y_help(ii).eq.0) krit=0 end do end do krit=0 j=1 do if(krit.ne.0)exit a: if (y(j).eq.0) then b: if (y(j).eq.y(j+1)) then jj=j+1 qwer=size(y)-1 allocate (vecout(qwer)) call schnitt(y,jj,vecout,qwer) deallocate (y) allocate (y(qwer)) y=vecout deallocate (vecout) jj=j+1 qwer=size(x)-1 allocate (vecout(qwer)) call schnitt(x,jj,vecout,qwer) deallocate (x) allocate (x(qwer)) x=vecout deallocate (vecout) x(j)=x(j)+1 else b j=j+1 endif b else a j=j+1 endif a if (j.ge.size(y)) krit=1 end do do j=1,size(y) if (y(j).eq.0) x(j)=x(j)+1 end do j=1 do if(j.ge.(size(y)-1))exit if (y(j).eq.0) then l=j+1 do l=l+1 if (l.ge.size(y)) exit if ((y(l).eq.0)) exit end do if (y(l).eq.0) then jj=j+1 qwer=size(y)-1 allocate (vecout(qwer)) call schnitt(y,jj,vecout,qwer) deallocate (y) allocate (y(qwer)) y=vecout deallocate (vecout) jj=j+1 qwer=size(x)-1 allocate (vecout(qwer)) call schnitt(x,jj,vecout,qwer) deallocate (x) allocate (x(qwer)) x=vecout deallocate (vecout) endif endif j=j+1 end do if (i.eq.row) then allocate (yy(size(y))) allocate (xx(size(x))) yy=y xx=x deallocate (y) deallocate (x) allocate (y(size(yy)+1)) allocate (x(size(xx)+1)) x=xx y=yy y(size(yy)+1)=1 x(size(xx)+1)=1 deallocate (yy,xx) end if allocate (c(s+1)) allocate (con(s+1)) c(1)=ssd(1)/n_gew(1) do j=2,s+1 nsum=sum(c(1:j-1)*n_gew(1:j-1)) c(j)=(ssd(j)-nsum)/n_gew(j) end do con=c mmm=x(1) cm(i,1:mmm)=con(1) do j=2,size(y) nsum=sum(x(1:j-1)) mm1=nsum+1 mm2=nsum+x(j) cm(i,mm1:mm2)=con(j) end do deallocate (ssd,nn,nn_help) deallocate (y_help,n_gew,ssd_help) deallocate (x,y) deallocate(c,con) end do C do ii=1,anz-1 C write(6,*)(cm(ii,j),j=1,k) C end do deallocate(mat2,mat,n1,n2) do j=1,k do i=1,mm cm(i,j)=cm(i,j)*znn(j) end do end do return end subroutine schnitt(vec,stelle,vecout2,qwer) C SCHNITT cuts the place STELLE out of the vector VEC C and returns a vector whose dimension is reduced by 1 integer stelle,qwer real,dimension(qwer+1)::vec real,dimension(qwer)::vecout2 if (stelle.eq.qwer+1) then vecout2=vec(1:qwer) else vecout2(1:stelle)=vec(1:stelle) vecout2(stelle:qwer)=vec(stelle+1:qwer+1) endif return end Subroutine dist212(dis1) C Subroutine calculates inverse distances from the origin to C boundaries of convex region (swtypes=21,22) use comm dis1=zero do i=1,mm a=tte(i)/step(irank(i)) if(a.gt.dis1)dis1=a end do return end SUBROUTINE indexx(n,arr,indx) INTEGER n,indx(n),M,NSTACK REAL arr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK) REAL a do 11 j=1,n indx(j)=j 11 continue jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 13 j=l+1,ir indxt=indx(j) a=arr(indxt) do 12 i=j-1,l,-1 if(arr(indx(i)).le.a)goto 2 indx(i+1)=indx(i) 12 continue i=l-1 2 indx(i+1)=indxt 13 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 itemp=indx(k) indx(k)=indx(l+1) indx(l+1)=itemp if(arr(indx(l)).gt.arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endif if(arr(indx(l+1)).gt.arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir) indx(ir)=itemp endif if(arr(indx(l)).gt.arr(indx(l+1)))then itemp=indx(l) indx(l)=indx(l+1) indx(l+1)=itemp endif i=l+1 j=ir indxt=indx(l+1) a=arr(indxt) 3 continue i=i+1 if(arr(indx(i)).lt.a)goto 3 4 continue j=j-1 if(arr(indx(j)).gt.a)goto 4 if(j.lt.i)goto 5 itemp=indx(i) indx(i)=indx(j) indx(j)=itemp goto 3 5 indx(l+1)=indx(j) indx(j)=indxt jstack=jstack+2 if(jstack.gt.NSTACK)then write(6,*)'NSTACK too small in indexx' write(8,*)'NSTACK too small in indexx' endif if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END SUBROUTINE rank(n,indx,irank) INTEGER n,indx(n),irank(n) INTEGER j do 11 j=1,n irank(indx(j))=j 11 continue return END