c This program performs sample size calculations
c for balanced ANOVAs
c in association with the power Web program.
c Currently, it will handle at most 5 factors, at most 5
c levels per factor, and at most 5000 observations.
c
c It will also handle at most 5 cb1's, 5 cv's, and 5 sigma's.
c Steve Verrill, 3/14/94
c Steve Verrill, 11/16/95
c Steve Verrill, 1/11/97
c The parameter and dimension statements in the anocpow, hyptest,
c power, and kofn subroutines must be in accord.
parameter (maxfac=5,maxlev=25,mxlevp1=maxlev+1,
x mxlevm1=maxlev-1,nrdx=5000,
x ncdx=mxlevp1,ncdc=mxlevm1)
real temp(ncdx)
real pred(nrdx),predsort(nrdx)
c Note that a copy of xmat is created in g in hyptest.
real*8 xmat(nrdx,ncdx)
real*8 y(nrdx)
real*8 c(ncdx,ncdc)
real*8 hypss
c??? real*8 resss ??? never assigned a value or used ???
real fracs(maxlev)
real*8 cb1(5),sigma(5),cb(mxlevm1),powerl
real*8 alpha
c
c If I were more careful,
c I could probably use fewer id vectors.
c
integer iorigord(nrdx)
character*20 id(nrdx),idorig(nrdx)
integer idoink(nrdx)
integer id2(nrdx)
character*20 idbase(nrdx)
integer idlevel(nrdx,maxfac)
integer numlev(maxfac)
character*80 datafile,results
character*10 level(maxfac,maxlev)
character*80 argv,intwrite
common /alphac/ alpha(3)
common /nrowc/ nrow
common /cvc/ cv(5)
call getarg(10,results)
open(unit=7,file=results,access="append",status="unknown")
call getarg(11,results)
write(6,9900) results
9900 format("The results from this run have been appended to",
x " the file ",a80,"
To see how to download",
x " this file,",
x ' read this description of anonymous ftp',
x '.','Alternatively, you should be able to do a',
x ' "Save As" or "Print" in your browser.
')
write(7,9700)
9700 format(////,1x,"**************************************",//,1x,
x "Output from power.2.f",//,1x,
x "**************************************",/)
c
c number of factors
c
call getarg(1,argv)
write(intwrite,7) argv
7 format(a80)
read(intwrite,*) numfac
c read(3,*) numfac
write(6,9501) numfac
9501 format("
The number of factors was ",i4,".
")
write(7,9601) numfac
9601 format(/,1x,'The number of factors was ',i2,'.',/)
call getarg(2,argv)
write(intwrite,7) argv
read(intwrite,*) (numlev(i),i = 1,numfac)
c read(3,*) (numlev(i),i = 1,numfac)
write(6,9502)
9502 format("The numbers of levels were:
")
write(7,9602) (numlev(i),i = 1,numfac)
9602 format(/,1x,'The numbers of levels were:',//,1x,5(i2,2x),/)
do 9504 i = 1,numfac
write(6,9503) numlev(i)
9503 format("
",i4)
9504 continue
call getarg(3,argv)
write(intwrite,7) argv
read(intwrite,*) numrep
c read(3,*) numrep
write(6,9505) numrep
9505 format("
The number of replicates was ",i4,
x".
")
write(7,9603) numrep
9603 format(/,1x,'The number of replicates was ',i4,'.',/)
c
c generate the 0,1 portion of the design matrix
c
c need dynamic memory allocation??? Yes --- probably don't need
c to handle more than 5,000 specimens, but should allow up to 10
c levels per factor (say) and at least five factors and 1000 reps
c which would require ...
if (numfac .eq. 1) then
nrow = 0
nlevtot = numlev(1)
do 150 i = 1,numlev(1)
do 142 j = 1,numlev(1)
temp(j) = 0.0
142 continue
temp(i) = 1.0
do 146 ii = 1,numrep
nrow = nrow + 1
do 144 j = 1,numlev(1)
xmat(nrow,j) = temp(j)
144 continue
146 continue
150 continue
elseif (numfac .eq. 2) then
nrow = 0
nlevtot = numlev(1) + numlev(2)
do 170 i = 1,numlev(1)
do 168 j = 1,numlev(2)
do 152 k = 1,nlevtot
temp(k) = 0.0
152 continue
temp(i) = 1.0
temp(numlev(1)+j) = 1.0
do 156 ii = 1,numrep
nrow = nrow + 1
do 154 k = 1,nlevtot
xmat(nrow,k) = temp(k)
154 continue
156 continue
168 continue
170 continue
elseif (numfac .eq. 3) then
c write(6,9506)
c 9506 format("
Got to 9506.
")
nrow = 0
nlevtot = numlev(1) + numlev(2) + numlev(3)
do 200 i = 1,numlev(1)
do 198 j = 1,numlev(2)
do 196 k = 1,numlev(3)
do 172 l = 1,nlevtot
temp(l) = 0.0
172 continue
temp(i) = 1.0
temp(numlev(1)+j) = 1.0
temp(numlev(1)+numlev(2)+k) = 1.0
do 176 ii = 1,numrep
nrow = nrow + 1
do 174 l = 1,nlevtot
xmat(nrow,l) = temp(l)
174 continue
176 continue
196 continue
198 continue
200 continue
elseif (numfac .eq. 4) then
nrow = 0
nlevtot = numlev(1) + numlev(2) + numlev(3) + numlev(4)
do 240 i = 1,numlev(1)
do 239 j = 1,numlev(2)
do 238 k = 1,numlev(3)
do 237 l = 1,numlev(4)
do 202 m = 1,nlevtot
temp(m) = 0.0
202 continue
temp(i) = 1.0
temp(j+numlev(1)) = 1.0
temp(k+numlev(1)+numlev(2)) = 1.0
temp(l+numlev(1)+numlev(2)+numlev(3)) = 1.0
do 206 ii = 1,numrep
nrow = nrow + 1
do 204 m = 1,nlevtot
xmat(nrow,m) = temp(m)
204 continue
206 continue
237 continue
238 continue
239 continue
240 continue
else
c
c numfac = 5
c
nrow = 0
nlevtot = numlev(1)+numlev(2)+numlev(3)+numlev(4)+numlev(5)
do 290 i = 1,numlev(1)
do 289 j = 1,numlev(2)
do 288 k = 1,numlev(3)
do 287 l = 1,numlev(4)
do 286 m = 1,numlev(5)
do 242 mm = 1,nlevtot
temp(mm) = 0.0
242 continue
temp(i) = 1.0
temp(j+numlev(1)) = 1.0
temp(k+numlev(1)+numlev(2)) = 1.0
temp(l+numlev(1)+numlev(2)+numlev(3)) = 1.0
temp(m+numlev(1)+numlev(2)+numlev(3)+numlev(4))
x = 1.0
do 246 ii = 1,numrep
nrow = nrow + 1
do 244 mm = 1,nlevtot
xmat(nrow,mm) = temp(mm)
244 continue
246 continue
286 continue
287 continue
288 continue
289 continue
290 continue
endif
c
c do the power calculations
c
ncol = nlevtot
3004 alpha(1) = .10
alpha(2) = .05
alpha(3) = .01
call getarg(4,argv)
write(intwrite,7) argv
read(intwrite,*) idfac
write(6,9506) idfac
9506 format("The factor tested was ",i1,".
",
x"This was a test for a difference in two",
x" pre-specified levels.
")
write(7,9604) idfac
9604 format(/,1x,'The factor tested was ',i1,'.',///,
x1x,'This was a test for a difference in two',
x' pre-specified levels.',/)
c read(3,*) idfac
c
c comparing 2 pre-specified levels --- idtest = 0
c
idtest = 0
c 3044 write(6,3045)
c write(7,3045)
c 3045 format(/,1x,'What are the 2 levels that you wish to compare?',
c x/)
c read(5,*) idlev1,idlev2
c write(7,*) idlev1,idlev2
c ****************************************************
c
c Since the design is balanced. The power
c should not depend on which 2 levels of the factor
c are chosen.
c
c*****************************************************
idlev1 = 1
idlev2 = 2
c if ((idlev1 .lt. 1) .or. (idlev1 .gt. numlev(idfac))) then
c write(6,3046) numlev(idfac)
c write(7,3046) numlev(idfac)
c 3046 format(/,1x,'The levels must lie between 1 and ',i2,/)
c go to 3044
c elseif ((idlev2 .lt. 1) .or. (idlev2 .gt. numlev(idfac))) then
c write(6,3046) numlev(idfac)
c write(7,3046) numlev(idfac)
c go to 3044
c elseif (idlev1 .eq. idlev2) then
c write(6,3047)
c write(7,3047)
c 3047 format(/,1x,'You must specify 2 distinct levels.',/)
c go to 3044
c endif
do 3048 i = 1,ncol
c(i,1) = 0.0
3048 continue
numpre = 0
do 3049 i = 1,idfac-1
numpre = numpre + numlev(i)
3049 continue
i1 = numpre + idlev1
i2 = numpre + idlev2
c(i1,1) = 1
c(i2,1) = -1
c 3052 write(6,3054)
c write(7,3054)
c 3054 format(/,1x,'How many differences do you want to consider?',
c x' (5 or fewer)',/)
c read(5,*) numdelta
c write(7,*) numdelta
call getarg(5,argv)
write(intwrite,7) argv
read(intwrite,*) numdelta
numcb1 = numdelta
c if ((numdelta .lt. 1) .or. (numdelta .gt. 5)) then
c write(6,3056)
c write(7,3056)
c 3056 format(/,1x,'The number must lie between 1 and 5',
c x ' (inclusive).',/)
c go to 3052
c endif
call getarg(6,argv)
write(intwrite,7) argv
read(intwrite,*) (cb1(i),i = 1,numdelta)
c write(6,3060)
c write(7,3060)
c 3060 format(/,1x,'What are they? (e.g., .20 for a 20% difference',
c x ' in means)',/)
c read(5,*) (cb1(i),i = 1,numdelta)
c write(7,*) (cb1(i),i = 1,numdelta)
nq = 1
call getarg(7,argv)
write(intwrite,7) argv
read(intwrite,*) numcv
c read(3,*) numcv
numsig = numcv
call getarg(8,argv)
write(intwrite,7) argv
read(intwrite,*) (cv(i),i = 1,numcv)
c read(3,*) (cv(i),i = 1,numcv)
call getarg(9,argv)
write(intwrite,7) argv
read(intwrite,*) powerl
write(7,9605) powerl
9605 format(/,1x,'The power desired was ',f5.3,'.',/)
c read(3,*) powerl
c
c prepare for the power calculations
c
call hyptest(xmat,y,nrdx,nrow,ncol,c,ncdx,nq,hypss,ihypdof,
1 notestf,2)
c the notestf (not estimable flag) should never be anything other than 0
c since we are only feeding it standard designs and considering standard
c tests, but could still test for notestf > 0
do 5810 jjj = 1,numcv
sigma(jjj) = cv(jjj)
5810 continue
call power(nq,cb1,numcb1,sigma,numsig,cb,powerl,idtest)
stop
end
c
c matrix addition
c
subroutine matadd(a,b,c,nrow,ncol,nrdim)
implicit double precision (a-h,o-z)
dimension a(nrdim,1),b(nrdim,1),c(nrdim,1)
do 20 i = 1,nrow
do 10 j = 1,ncol
c(i,j) = a(i,j) + b(i,j)
10 continue
20 continue
return
end
c
c matrix subtraction
c
subroutine matsub(a,b,c,nrow,ncol,nrdim)
implicit double precision (a-h,o-z)
dimension a(nrdim,1),b(nrdim,1),c(nrdim,1)
do 20 i = 1,nrow
do 10 j = 1,ncol
c(i,j) = a(i,j) - b(i,j)
10 continue
20 continue
return
end
c
c matrix multiplication
c
subroutine matmult(a,nra,nca,b,nrb,ncb,c,nrda,nrdb)
implicit double precision (a-h,o-z)
dimension a(nrda,1),b(nrdb,1),c(nrda,1)
if (nca .ne. nrb) then
write(6,10)
c write(7,10)
10 format(1x,'Matrix multiplication error:',/,
1 1x,'The a column and b row dimensions do not match.')
stop
else
do 50 j = 1,ncb
do 40 i = 1,nra
c(i,j) = 0.0
do 30 k = 1,nca
c(i,j) = c(i,j) + a(i,k)*b(k,j)
30 continue
40 continue
50 continue
endif
return
end
c
c matrix transposition
c
subroutine mattran(a,nra,nca,at,nrda,ncda)
implicit double precision (a-h,o-z)
dimension a(nrda,1),at(ncda,1)
do 20 i = 1,nra
do 10 j = 1,nca
at(j,i) = a(i,j)
10 continue
20 continue
return
end
function dot(x,y,n)
implicit double precision (a-h,o-z)
dimension x(n),y(n)
dot = 0.0
do 10 i = 1,n
dot = dot + x(i)*y(i)
10 continue
return
end
function vnorm(x,n)
implicit double precision (a-h,o-z)
dimension x(n)
external dot
vnorm = sqrt(dot(x,x,n))
return
end
subroutine vector(c,nrow,j,v,nrdc)
c This places the jth column of the matrix, c,
c into vector v.
implicit double precision (a-h,o-z)
dimension c(nrdc,1),v(nrow)
do 10 i = 1,nrow
v(i) = c(i,j)
10 continue
return
end
subroutine hyptest(x,y,nrdx,nrow,ncol,c,ncdx,nq,hypss,ihypdof,
1 notestf,idesign)
c If idesign equals 0,
c this subroutine checks whether the c's are estimable and
c if they are it tests the hypothesis c1'b = ... = cq'b = 0.
c If idesign equals 1 (estimability calculations requested) or
c idesign equals 2 (power calculations requested),
c this subroutine checks the estimability of the c's and provides for
c power calculations.
c The parameter and dimension statements in the anocpow, hyptest,
c power, and kofn subroutines must be in accord.
parameter (maxfac=5,maxlev=25,
x nrowd=5000,
x nrdc=maxlev+1)
implicit double precision (a-h,o-z)
dimension x(nrdx,ncol),y(nrdx),c(ncdx,ncol),e(nrdc),work(nrowd)
dimension f(nrdc,nrdc),ft(nrdc,nrdc),fft(nrdc,nrdc)
dimension vtemp(nrdc),temp1(nrdc),temp2(nrdc)
dimension g(nrowd,nrdc),h(nrowd,nrdc)
dimension ua(nrowd,nrdc),xlama(nrdc),va(nrdc,nrdc)
c What is notest? It stands for "not estimable". It is
c a listing of contrasts that do not appear to be estimable.
c notestf lets us know whether there any non-estimable contrasts.
dimension notest(nrdc)
dimension at(nrdc,nrowd),ata(nrdc,nrdc)
dimension det(2)
c****** this common statement is unneccessary? need to dimension
c****** but not in a common statement????
common /hec/ ux(nrowd,nrdc),vx(nrdc,nrdc),xlamx(nrdc),
x a(nrowd,nrdc)
common /ataic/ atai(nrdc,nrdc)
common /ixrankc/ ixrank
external vnorm
c*** I'm almost certainly be wasteful of memory in this and
c*** associated routines, but I don't currently have the
c*** time or the incentive to clean it up. It works fine
c*** on the machines that I use. 5000x26 --- lots of these
c*** so why aren't I running out of space????? Would
c*** on a PC?
c determining machine double precision
xmult = .99d0
delta = 1.0d0
2 delta = xmult*delta
if ((1.0d0 + delta) .gt. 1.0d0) go to 2
delta = delta/xmult
threshm = sqrt(delta)
notestf = 0
do 5 i = 1,nrow
do 4 j = 1,ncol
g(i,j) = x(i,j)
4 continue
c write(7,*) (g(i,j),j = 1,ncol)
5 continue
c
c ssvdc is the linpack singular value decomposition routine
c
call dsvdc(g,nrdx,nrow,ncol,xlamx,e,ux,nrdx,vx,ncdx,
xwork,21,info)
c
c If info is not equal to 0 print an error message and stop.
c
if (info .ne. 0) then
write(6,7) info
c write(7,7) info
7 format('
ERROR 7: The info value from a singular',
x ' value decomposition was ',i3,'.
Call Steve',
x ' Verrill at 608-231-9375 for assistance.
')
stop
endif
c
c obtain the threshold that is used to determine whether singular values
c are "essentially zero"
c
xnorm = 0.0d0
do 10 i = 1,ncol
xnorm = xnorm + xlamx(i)*xlamx(i)
10 continue
xnorm = sqrt(xnorm)
thresh = xnorm*threshm
c
c determine the rank of x
c
do 120 i = 1,ncol
if (xlamx(i) .gt. thresh) go to 120
ixrank = i - 1
go to 130
120 continue
ixrank = ncol
c
c determine whether the c's are estimable (there exists an a
c such that x'a = c)
c
130 iresdof = nrow - ixrank
do 150 i = 1,ncol
do 140 j = 1,ixrank
f(i,j) = vx(i,j)
140 continue
150 continue
c
c obtain the transpose of f
c
c******* could just call mattran
do 22 i = 1,ncol
do 12 j = 1,ixrank
ft(j,i) = f(i,j)
12 continue
22 continue
c
c multiply f by ft to obtain fft
c
c******* could just call matmult
do 52 j = 1,ncol
do 42 i = 1,ncol
fft(i,j) = 0.0d0
do 32 k = 1,ixrank
fft(i,j) = fft(i,j) + f(i,k)*ft(k,j)
32 continue
42 continue
52 continue
c
c fft is the projection operator onto the space spanned by
c the columns of v (or the columns of x')
c
mm = 0
do 200 i = 1,nq
c******* could replace vnorm with BLAS dnrm2
call vector(c,ncol,i,vtemp,ncdx)
thresh = vnorm(vtemp,ncol)*threshm
call matmult(fft,ncol,ncol,vtemp,ncol,1,temp1,ncdx,ncdx)
call matsub(vtemp,temp1,temp2,ncol,1,ncdx)
if (vnorm(temp2,ncol) .lt. thresh) go to 200
mm = mm + 1
notest(mm) = i
200 continue
if (mm .eq. 0) go to 300
c****** note that the stuff below only allows 99 c's (i2 write format)
write(6,210) (notest(i),i = 1,mm)
write(7,9210) (notest(i),i = 1,mm)
210 format("
The following c's do not appear to be",
x" estimable",
1"
",(1x,10(i2,2x)),"
")
9210 format(///,1x,"The following c's do not appear to be estimable",
1//,(1x,10(i2,2x)),/)
notestf = 1
return
c
c calculate the ai's that satisfy x'ai = ci
c
300 if (idesign .ne. 1) go to 305
write(6,302)
write(7,9302)
302 format("
The c's all appear to be estimable.
")
9302 format(///,1x,"The c's all appear to be estimable.",/)
305 do 320 i = 1,nrow
do 310 j = 1,ixrank
g(i,j) = ux(i,j)/xlamx(j)
310 continue
320 continue
c
c h is the matrix that transforms c's into a's
c (h = u(nxr)*diaginv(rxr)*vt(rxp))
c
call matmult(g,nrow,ixrank,ft,ixrank,ncol,h,nrdx,ncdx)
call matmult(h,nrow,ncol,c,ncol,nq,a,nrdx,ncdx)
if (idesign .eq. 0) go to 323
call mattran(a,nrow,nq,at,nrdx,ncdx)
call matmult(at,nq,nrow,a,nrow,nq,ata,ncdx,nrdx)
c
c Taken together, dpofa and dpodi yield an inverse of the
c ata matrix. dpofa makes use of the Cholesky algorithm
c to factor the positive definite matrix ata.
c dpofa and dpodi are linpack routines. Their source is
c included in this distribution. See the linpack manual
c for details. (Dongarra, Moler, Bunch, and Stewart,
c Linpack Users' Guide, SIAM, 1979)
c
call dpofa(ata,ncdx,nq,info)
if (info .ne. 0) then
write(6,9999) info
c write(7,9999) info
9999 format('
ERROR 9999: The info value from a matrix',
x ' inversion was',1x,i3,'
Call Steve Verrill at',
x ' 608-231-9375 for assistance.
')
stop
endif
c
c calculate the inverse of ata
c
job = 1
call dpodi(ata,ncdx,nq,det,job)
do 322 i = 1,nq
do 321 j = i,nq
atai(i,j) = ata(i,j)
atai(j,i) = ata(i,j)
321 continue
322 continue
return
c
c idesign = 0
c
323 do 325 i = 1,nrow
do 324 j = 1,nq
g(i,j) = a(i,j)
324 continue
325 continue
c
c print the a vectors
c
c write(8,326)
c326 format(//,1x,"The a's are (as columns)",/)
c do 328 i = 1,nrow
c write(8,327) (a(i,j),j = 1,nq)
c327 format(10(1x,e11.3))
c328 continue
c
c perform a singular value decomposition of a
c in order to obtain an orthonormal basis of L(a)
c (the linear span of the columns of a)
c
call dsvdc(g,nrdx,nrow,nq,xlama,e,ua,nrdx,va,ncdx,
xwork,20,infoa)
c It appears to me that can get the non-centrality parameter without
c dpofa and dpodi if have ua --- no if only know c'b values
c but yes if we know Xb values.
if (infoa .ne. 0) then
write(6,327) infoa
c write(7,327) infoa
327 format('
ERROR 327: The info value from a singular',
x ' value decomposition was ',i3,'.
Call Steve',
x ' Verrill at 608-231-9375 for assistance.
')
stop
endif
c
c obtain the hypothesis sum of squares
c
hypss = 0.0d0
do 340 j = 1,nq
temp0 = 0.0d0
do 330 i = 1,nrow
temp0 = temp0 + ua(i,j)*y(i)
330 continue
hypss = hypss + temp0*temp0
340 continue
ihypdof = nq
c
c obtain the residual sum of squares
c
xss = 0.0d0
do 440 j = 1,ixrank
temp0 = 0.0d0
do 430 i = 1,nrow
temp0 = temp0 + ux(i,j)*y(i)
430 continue
xss = xss + temp0*temp0
440 continue
totss = 0.0d0
do 450 i = 1,nrow
totss = totss + y(i)*y(i)
450 continue
c 11/18/95 resss never returned ?
resss = totss - xss
return
end
c
c
c
subroutine power(nq,cb1,numcb1,sigma,numsig,cb,powerl,idtest)
c
c This subroutine performs power calculations.
c
c The parameter and dimension statements in the anocpow, hyptest,
c power, and kofn subroutines must be in accord.
parameter (maxlev=25,mxlevm1=maxlev-1,ncdx=maxlev+1)
implicit double precision (a-h,o-z)
dimension cb1(5),sigma(5),cb(mxlevm1)
c What is the 5,5,3?
c first five -- sigmas second five -- cb1's three -- alphas
dimension powera(5,5,3)
dimension numrep(5,5,3)
c What is the 3,5? 3 -- alphas 5 -- sigmas
dimension halfwidt(3,5)
real cv
common /alphac/ alpha(3)
common /ataic/ atai(ncdx,ncdx)
common /nrowc/ nrow
common /ixrankc/ ixrank
common /cvc/ cv(5)
iresdof = nrow - ixrank
if (idtest .eq. 1) go to 1000
c
c comparison of two particular levels of a factor
c
d = sqrt(atai(1,1))
do 200 ii = 1,3
f = iresdof
c
c cdft is a routine from DCDFLIB. It performs calculations
c associated with the t distribution. DCDFLIB is a
c public domain library of "routines for cumulative distribution
c functions, their inverses, and their parameters." It was produced
c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and
c Kathy Russell of the Department of Biomathematics,
c M.D. Anderson Cancer Center, The University of Texas.
c One way to obtain it is through statlib.
c To get started on statlib, send the e-mail message
c
c send index
c
c to statlib@lib.stat.cmu.edu.
c The statlib URL is
c http://lib.stat.cmu.edu/
c or, if that doesn't work
c http://lib.stat.cmu.edu/.index.html
c
p = 1.0 - alpha(ii)/2.0
call cdft(2,p,1.0d0-p,xbase,f,istatus,bound)
if (istatus .ne. 0) then
write(6,10) istatus
c write(7,10) istatus
10 format('
The istatus value from a T inverse',
x ' routine was ',i2,'.
Please call Steve Verrill',
x ' at 608-231-9375 for assistance.
')
stop
endif
do 100 i = 1,numsig
do 90 j = 1,numcb1
xncp = (cb1(j)*d)/sigma(i)
c
c cdffnc is a routine from DCDFLIB. It performs calculations
c associated with the non-central F distribution. DCDFLIB is a
c public domain library of "routines for cumulative distribution
c functions, their inverses, and their parameters." It was produced
c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and
c Kathy Russell of the Department of Biomathematics,
c M.D. Anderson Cancer Center, The University of Texas.
c One way to obtain it is through statlib.
c To get started on statlib, send the e-mail message
c
c send index
c
c to statlib@lib.stat.cmu.edu.
c The statlib URL is
c http://lib.stat.cmu.edu/
c or, if that doesn't work
c http://lib.stat.cmu.edu/.index.html
c
f = iresdof
call cdffnc(1,p,q,xbase*xbase,1.0d0,f,xncp*xncp,
x istatus,bound)
if (istatus .ne. 0) then
write(6,1025) istatus
c write(7,1025) istatus
stop
endif
powera(i,j,ii) = 1.0d0 - p
c
c find the number of replications needed to achieve the
c desired power ("1 replication" means the original design,
c "k replications" means k copies of the original design)
c
numrep(i,j,ii) = 1
if (powerl .le. powera(i,j,ii)) go to 90
numl = 1
numh = 50
xncpt = xncp*sqrt(50.0d0)
idof = 50*nrow - ixrank
f = idof
c
c cdft is a routine from DCDFLIB. It performs calculations
c associated with the t distribution. DCDFLIB is a
c public domain library of "routines for cumulative distribution
c functions, their inverses, and their parameters." It was produced
c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and
c Kathy Russell of the Department of Biomathematics,
c M.D. Anderson Cancer Center, The University of Texas.
c One way to obtain it is through statlib.
c To get started on statlib, send the e-mail message
c
c send index
c
c to statlib@lib.stat.cmu.edu.
c The statlib URL is
c http://lib.stat.cmu.edu/
c or, if that doesn't work
c http://lib.stat.cmu.edu/.index.html
c
p = 1.0 - alpha(ii)/2.0
call cdft(2,p,1.0d0-p,x,f,istatus,bound)
if (istatus .ne. 0) then
write(6,10) istatus
c write(7,10) istatus
stop
endif
c
c cdffnc is a routine from DCDFLIB. It performs calculations
c associated with the non-central F distribution. DCDFLIB is a
c public domain library of "routines for cumulative distribution
c functions, their inverses, and their parameters." It was produced
c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and
c Kathy Russell of the Department of Biomathematics,
c M.D. Anderson Cancer Center, The University of Texas.
c One way to obtain it is through statlib.
c To get started on statlib, send the e-mail message
c
c send index
c
c to statlib@lib.stat.cmu.edu.
c The statlib URL is
c http://lib.stat.cmu.edu/
c or, if that doesn't work
c http://lib.stat.cmu.edu/.index.html
c
call cdffnc(1,p,q,x*x,1.0d0,f,xncpt*xncpt,
x istatus,bound)
if (istatus .ne. 0) then
write(6,1025) istatus
c write(7,1025) istatus
stop
endif
pt = 1.0d0 - p
if (pt .ge. powerl) go to 50
c
c numrep is set to 100 so that ** is printed in the power table.
c ** indicates that the necessary number of replications is
c greater than 50.
c
numrep(i,j,ii) = 100
go to 90
50 idel = numh - numl
if (idel .lt. 2) then
numrep(i,j,ii) = numh
go to 90
endif
numt = numl + idel/2
xnumt = numt
xncpt = xncp*sqrt(xnumt)
idof = numt*nrow - ixrank
f = idof
c
c cdft is a routine from DCDFLIB. It performs calculations
c associated with the t distribution. DCDFLIB is a
c public domain library of "routines for cumulative distribution
c functions, their inverses, and their parameters." It was produced
c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and
c Kathy Russell of the Department of Biomathematics,
c M.D. Anderson Cancer Center, The University of Texas.
c One way to obtain it is through statlib.
c To get started on statlib, send the e-mail message
c
c send index
c
c to statlib@lib.stat.cmu.edu.
c The statlib URL is
c http://lib.stat.cmu.edu/
c or, if that doesn't work
c http://lib.stat.cmu.edu/.index.html
c
p = 1.0 - alpha(ii)/2.0
call cdft(2,p,1.0d0-p,x,f,istatus,bound)
if (istatus .ne. 0) then
write(6,10) istatus
c write(7,10) istatus
stop
endif
c
c cdffnc is a routine from DCDFLIB. It performs calculations
c associated with the non-central F distribution. DCDFLIB is a
c public domain library of "routines for cumulative distribution
c functions, their inverses, and their parameters." It was produced
c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and
c Kathy Russell of the Department of Biomathematics,
c M.D. Anderson Cancer Center, The University of Texas.
c One way to obtain it is through statlib.
c To get started on statlib, send the e-mail message
c
c send index
c
c to statlib@lib.stat.cmu.edu.
c The statlib URL is
c http://lib.stat.cmu.edu/
c or, if that doesn't work
c http://lib.stat.cmu.edu/.index.html
c
call cdffnc(1,p,q,x*x,1.0d0,f,xncpt*xncpt,
x istatus,bound)
if (istatus .ne. 0) then
write(6,1025) istatus
c write(7,1025) istatus
stop
endif
pt = 1.0d0 - p
if (pt .ge. powerl) go to 70
numl = numt
go to 50
70 numh = numt
go to 50
90 continue
c
c calculate approximate half-widths of confidence intervals
c
halfwidt(ii,i) = xbase*sigma(i)/d
100 continue
200 continue
write(6,205)
write(7,9205)
205 format('
Power tables
')
9205 format(///,1x,'Power tables',//)
do 300 ii = 1,3
write(6,310) alpha(ii)
write(7,9310) alpha(ii)
310 format('
The significance level is ',f4.2,
x '
')
9310 format(///,1x,'The significance level is ',f4.2,//)
write(6,320)
write(7,9320)
320 format('
The differences are on top.
',
x 'The coefficients of variation are along the left side.',
x '
')
9320 format(/,1x,'The differences are on top.',/,1x,
x 'The coefficients of variation are along the left side.',//)
write(6,8010)
8010 format("| ") do 330 k = 1,numcb1 write(6,8014) cb1(k) 8014 format(" | ",f5.3) 330 continue write(7,9330) (cb1(i),i = 1,numcb1) 9330 format(/,18x,5(e12.4,3x)) write(7,9332) 9332 format(/) do 340 i = 1,numsig write(6,8020) cv(i) 8020 format(" |
|---|---|
| ",f5.3) do 335 j = 1,numcb1 write(6,8030) powera(i,j,ii) 8030 format(" | ",f5.3) 335 continue 340 continue do 9340 i = 1,numsig write(7,9335) cv(i),(powera(i,j,ii),j = 1,numcb1) 9335 format(1x,e12.4,5x,5(e12.4,3x)) 9340 continue write(6,8040) 8040 format(" |