C MVI3.FOR Revised April 9, 2001 C Only change is name for Marsaglia's gamma generator. C Name changed from gamma to gamgen (lines 933,1190,1199). C Last revised October 8,1998. C MVI3.FOR is a Fortran 90 program to evaluate multivariate normal and multi- C variate-t integrals over a convex region bounded by hyperplanes. C It is assumed that the variance covariance (or correlation) matrix is C known except for a scalar constant (sigmasq). C If the scalar constant is known, the integral is multivariate normal. If C the scalar constant is unknown but the estimate has the usual chi-squared C relationship with associated degrees of freedom, then the integral is the C multivariate-t. C C The program requires that two files exist prior to running the program. C They are qcalc.in, which contains all input requirements, and qcalc.out C which will contain the output generated by the program. C C Contents of input file qcalc.in: C Line 1 contains 6 integers in free format: k, iseed, ndenom, mocar, irep, C mm. C The items are (in proper order): C k dimension of integral C iseed Integer in range 1 to 2**31 -1 for a 32 bit machine C ndenom degrees of freedom for variance estimate (-1 if var known) C mocar # of random directions used for each individual estimate C irep # of individual estimates C mm # of hyperplane boundaries for convex region C C Lines 2 to k + 1 contain the lower triangular portion of the variance C covariance (or correlation) matrix. Each line contains one row of the C matrix. C Lines k+2 to k+mm+1 contain the mm inequalities corresponding to the C convex region Lx.le.d where L is (mm, k), and x and d are vectors. C The ith inequality is represented as L(i,1) ... L(i,k) d(i). C C Several problems may be input sequentially, each problem requiring C k+mm+1 lines. Values for k and mm do not need to be the same for each C problem. C C The program expects a row of 4 or more -1's after the last problem. C C The convex region need not contain the origin, although if the origin is C on a boundary, or the origin is exterior to the region, "crude Monte C Carlo" is always used for the evaluation. C If the origin is on a boundary, a preferred procedure is to replace the C value 0. of d with 0.0000001. C C An alternate mode (binning) which is faster, usually by a factor of C approximately 2, can be invoked by replacing the constant mm (the number C of hyperplanes boundaries) by -mm. The program uses nq bins. The number C of bins required depends on the shape of the region of integration. More C specifically, it is a function of the minimum distance of the boundary C from the origin. The program uses an empirical formula to determine the C appropriate value for the value of nq. If the minimum distance is C greater than .2, nq is set to 25. if the minimum distance is less than C .01, the program uses the non-binning mode. For distances between .01 C and .20, the formula gives the appropriate value for nq. The user who C desires to use a fixed value for nq needs only to replace the statement C "nqfix=1" with the two statements "nqfix=a" and "nq=b" where a is any C integer other than 1, and b is the desired value for nq. C The binning method has been shown to be especially useful in calculating C the percentage points of certain statistics. C The program may also be used to estimate the integral using "crude" Monte C Carlo. Replacing k, the dimension of the integral by -k causes the C program to use "crude" Monte Carlo. C EXAMPLE 1: C Suppose we have a 3-dimensional multivariate normal integral with C covariance matrix the identity matrix. We have ranges of -1 to 2 C for the first dimension, -2.1 to 1.4 for the second dimension and C -.5 to infinity for the third. C C QCALC.IN might contain the following 10 lines. C 3 457 -1 10000 1 5 C 1 C 0 1 C 0 0 1 C -1 0 0 1 C 1 0 0 2 C 0 -1 0 2.1 C 0 1 0 1.4 C 0 0 -1 .5 C -1 -1 -1 -1 -1 -1 C C The above input would use 10000 x 1 = 10000 random directions to C produce a single estimate for the integral. C We might replace the first line with C 3 457 -1 1000 10 5 C or C 3 457 -1 100 100 5 C The first replacement would average 10 estimates, each using 1000 C random directions, and the second replacement would average 100 C estimates, each using 100 random directions. Each replacement C would enable an estimate of the accuracy of the estimate for the C value of the integral. The original first line would provide no error C estimate. C The program output for the three different inputs follows: C (all examples were run using a Pentium 90 processor.) C date and time 19981008102336.830 C elapsed time in seconds is 2.000 C no of pops is 3 seed is 457 C df for var est is 30 C no of ran dir. is 10000 C value of integral is calculated 1 times C mean value is 5.003167E-01 C -------------------- C date and time 19981008102335.520 C elapsed time in seconds is 1.000 C no of pops is 3 seed is 457 C df for var est is 30 C no of ran dir. is 1000 C value of integral is calculated 10 times C mean value is 5.003167E-01 C standard error of mean is 5.807164E-04 C -------------------- C date and time 19981008102338.100 C elapsed time in seconds is .000 C no of pops is 3 seed is 457 C df for var est is 30 C no of ran dir. is 100 C value of integral is calculated 100 times C mean value is 5.003167E-01 C standard error of mean is 9.492344E-04 C -------------------- C Suppose we change the input so that there are 17 df and we require C 100 estimates each with 100 directions. C Line 1 of the input becomes C 3 457 17 100 100 5 C The output becomes: C The output using 100 estimates each with 100 directions is: C date and time 19971223171318.080 C elapsed time in seconds is .000 C no of pops is 3 seed is 457 C df for var est is 17 C no of ran dir. is 100 C value of integral is calculated 100 times C mean value is 4.928888E-01 C standard error of mean is 1.035211E-03 C -------------------- C For all of the above cases, the region contained the origin C The program uses "crude" Monte Carlo when the region does C not contain the origin, and in that case a standard error of C the mean is given even when there is only a "single estimate." C The question arises as to how many random directions are needed. C Of course, the larger the number of random directions, the smaller C will be the standard error of the mean. For almost all cases, C 10000 will be sufficient to give an se of .005. This includes C values of k as large as 20. The shape of the region affects the error. C The more the region resembles an hypersphere, (center at the origin), C the smaller the standard error will be. C Three additional examples are given below: C SAMPLE INPUT (in file qcalc.in) C 4 457 -1 100 100 5 C 1 C .5 1 C .5 .5 1 C .5 .5 .5 1 C 2 -1 0 0 1 C 1 0 -1 0 1 C 0 0 -1 1 1 C -1 -1 2 0 1 C -1 -1 -4 0 1 C 3 457 30 10000 1 -6 C 1 C 0 1 C 0 0 1 C .2865 -.2865 0 1 C .2865 0 -.2865 1 C 0 .2865 -.2865 1 C -.2865 .2865 0 1 C -.2865 0 .2865 1 C 0 -.2865 .2865 1 C -1 -1 -1 -1 -1 -1 C C 3 457 30 10000 1 6 C 1 C 0 1 C 0 0 1 C .2865 -.2865 0 1 C .2865 0 -.2865 1 C 0 .2865 -.2865 1 C -.2865 .2865 0 1 C -.2865 0 .2865 1 C 0 -.2865 .2865 1 C -1 -1 -1 -1 -1 -1 C Running the program should produce the following output: C C SAMPLE OUTPUT (in file qcalc.out) C date and time 19981008102338.700 C elapsed time in seconds is 1.000 C no of pops is 4 seed is 457 C variance assumed known C no of ran dir. is 100 C value of integral is calculated 100 times C mean value is 1.811718E-01 C standard error of mean is 1.828163E-03 C -------------------- C date and time 19981008102339.200 C Binning method used C elapsed time in seconds is .000 C no of pops is 3 seed is 457 C df for var est is 30 C no of ran dir. is 10000 C value of integral is calculated 1 times C mean value is 9.499553E-01 C -------------------- C date and time 19981008102339.420 C Crude Monte Carlo is used C elapsed time in seconds is 1.000 C no of pops is 3 seed is 457 C df for var est is 30 C no of random directions is 10000 C points = 20000.000000 C estimate of integral is 9.491000E-01 C standard error is 1.554175E-03 C -------------------- C The methodology of the program is described in Somerville "Numerical C Computation of Multivariate Normal and Mutivariate t Probabilities Over C Convex Regions", Journal of Computational and Graphical Statistics (to C appear 1998) C The present program has been extended to include convex regions which do not include C the origin and, for the "binning" method uses an empirical formula to C determine the appropriate value for nq, the number of bins. C A number of subroutines used in the program have been adapted from C Numerical Methods, by Press et al (Cambridge University Press). module comm character(len=28)elapse character(len=8)date character(len=10)time real,dimension(:,:),allocatable::vcc,cmc,gd,hm,bound,ttemp real,dimension(:,:),allocatable::zerobnd real,dimension(:),allocatable::dc,wl,pl,prop real, dimension(:),allocatable::qans real,dimension(:),allocatable::dmd integer,dimension(:),allocatable::ibin,ifacee,ifaceb real sb,zmd2,zmnd2,zmdzn,znd2,fconst,zm,zn,cnlogn integer k,ndenom,mocar,irep,mm,istart,iend,negs,mmmvib real zero,zone,two,pi,zoneplus,sx,ssx,zrep,sd,zmean real dend integer iran(97),iy,irr,jseed,ndzero,moncarl,nz,nw,nran real cc integer ncount1,ncount2,count_rate,count,ier,count_max,iseed integer nq parameter (mr=714025,ia=1366,ic=150889,rm=1./mr) end module comm use comm C uses portable random number generator C vc input var-covar matrix C qans vector of q values (irep of them) C ifaceb(i),ifacee(i) designation of ith pair of hyperplanes real se,seconds,points,poff,xxx elapse=' elapsed time in seconds is ' Open(8,file='qcalc.out',access='sequential',form='formatted', $recl=80,status='old') Open(9,file='qcalc.in',access='sequential',form='formatted', $recl=80,status='old') 21 read(9,*)k,iseed,ndenom,mocar,irep,mm write(6,*)k,iseed,ndenom,mocar,irep,mm moncarl=0 if((k.lt.0).and.(mocar.lt.0).and.(ndenom.lt.0).and. &(mocar.lt.0))go to 997 C If moncarl=1, crude Monte Carlo is used regardless of sign of mm. if(k.le.-1)then moncarl=1 k=abs(k) end if C following two statements were added so the program could use both C mvi and mvib methods to evaluate the integrals. mmmvib=mm mm=abs(mm) C set miscellaneous constants zero=0. zone=1. two=2. pi=3.14159265359 iseed=abs(iseed) jseed =iseed C iy set to a fixed value for random # generator iy=0 C fill random number bin do j=1,97 iseed=mod(ia*iseed+ic,mr) iran(j)=iseed enddo nz=iran(27) nw=iran(57) nran=iran(97) C irr is seed for racos irr=iseed zm=k 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(cmc(k,k) ) allocate( gd(mm,k)) allocate(dc(k),dmd(mm)) allocate(vcc(k,k)) C zero storage for q calculations allocate(qans(irep)) qans=0. C Calculate constant for recdist distribution if(ndenom.gt.0) fconst=two*beta(zmd2,znd2) C Input variance covariance matrix (lower triangular portion only) C Input one row at a time do i=1,k read(9,*)(vcc(i,j),j=1,i) end do do i=1,k-1 do j=i+1,k vcc(i,j)=vcc(j,i) end do end do C***************************************************** call chol C***************************************************** C write(8,*)'cholesky decomp' C do 7651 i=1,k C write(8,*)(cmc(i,j),j=1,k) C7651 continue C************************************************ deallocate(vcc) C Input boundaries in form Lx.le.d where the elements of d may be C positive, negative or zero. The program places the boundaries with C negative RHS elemnts in the first "negs" positions. The boundaries C with zero RHS are placed in positions negs+1 to negs+ ndzero. C Boundaries with +ive RHS are placed in the last mm-(negs+ndzero) C positions. For the boundaries with +ive or -ive RHS, the LHS is C scaled so that the RHS is +1 or -1 (sign implied by position). allocate(bound(mm,k),zerobnd(mm,k)) negs=0 ndzero=0 istart = 1 iend = mm+1 do i=1,mm read(9,*)(gd(i,j),j=1,k),dend if(dend.eq.zero) then ndzero=ndzero+1 zerobnd(ndzero,1:k) = gd(i,1:k) go to 110 endif if (dend.gt.zero)then iend =iend-1 bound(iend,1:k)=gd(i,1:k)/dend else dend=-dend negs=negs+1 bound(istart,1:k)=gd(i,1:k)/dend istart = istart+1 end if 110 end do if(ndzero.gt.0)then write(8,*)'Warning:RHS of',ndzero,' eqns eq zero' bound(negs+1:negs+ndzero,1:k)=zerobnd(1:ndzero,1:k) end if deallocate(zerobnd) C Transform boundaries to new coordinate system bound=matmul(bound,cmc) C write(8,*)'boundaries in new coord system' C do i=1,mm C write(8,*)(bound(i,j),j=1,k) C end do C normalize boundary eqns allocate (ttemp(mm,mm)) ttemp=matmul(bound,transpose(bound)) do i=1,mm ttemp(i,i)=sqrt(ttemp(i,i)) dmd(i)=zone/ttemp(i,i) do j=1,k gd(i,j)=bound(i,j)/ttemp(i,i) end do end do deallocate(ttemp) C dmd(i) is distance from the origin (Unless plane goes through origin) C Prepare for a sort which will order planes from closest to origin to C farthest from origin. Note: changing the sign of dmd(i) does not change C the pair of values returned by mvidist! do i=1,negs dmd(i)=-dmd(i) end do do i=negs+1,negs+ndzero dmd(i)=zero end do if(mm.gt.1)call sortary(mm,k,dmd,gd) call date_and_time(date,time) call system_clock(count,count_rate,count_max) ncount1=count write(8,*)'date and time ',date,time C If origin is not interior to the convex region or C Monte Carlo is the method of choice, use subroutine crude. if((moncarl.eq.1).or.(ndzero.gt.0).or.(negs.gt.0))then write(8,*)'Crude Monte Carlo is used' mocar=mocar*irep call crude go to 950 end if C If binning (MVIB) is the method of choice and mocar=0, and C origin is interior to convex region, then use binning method. if((moncarl.eq.0).and.(mmmvib.lt.0).and.(ndzero.eq.0).and. &(negs.eq.0))then call binning C Check to see if binning abandoned if(mmmvib.gt.0)then write(8,*)'Boundary too close to origin, binning not used' go to 580 end if go to 901 end if 580 do 900 njj=1,irep sum=0. if(ndenom.gt.0) then do 600 i=1,mocar call racos C find reciprocal distance to nearest face/mmmvib call mvidist(dis1,dis2) if(dis1.eq.0.) then sum=sum+zone go to 590 end if qsphere=zm*dis1*dis1 qsphere=zone/qsphere call pf(k,ndenom,qsphere,poff,ier) sum=sum + poff 590 if(dis2.eq.0.) then sum=sum+zone go to 600 end if qsphere=zm*dis2*dis2 qsphere=zone/qsphere call pf(k,ndenom,qsphere,poff,ier) sum=sum + poff 600 continue end if if(ndenom.lt.0) then xxx=zm/two do i=1,mocar call racos C find reciprocal distance to nearest face call mvidist(dis1,dis2) if(dis1.eq.0) then sum=sum+zone go to 605 end if qsphere=dis1*dis1*two qsphere=1/qsphere poff=gammp(xxx,qsphere) sum=sum+poff 605 if(dis2.eq.0) then sum=sum+zone go to 610 end if qsphere=dis2*dis2*two qsphere=1/qsphere poff=gammp(xxx,qsphere) sum=sum+poff 610 end do end if points=mocar*2 qans(njj)=sum/points cc=qans(1) cm=qans(njj)-cc sx=sx+cm ssx=ssx+cm*cm 900 continue 901 call date_and_time(date,time) call system_clock(count) ncount2=count seconds=ncount2-ncount1 if(seconds.lt.0)seconds=seconds+count_max seconds=seconds/count_rate write(8,8989)elapse,seconds 8989 format(a28,f10.3) sd=9999. se=9999. if(zrep.eq.1) go to 902 sd=(ssx-sx*sx/zrep)/(zrep-1.) se=sqrt(sd/zrep) sd=sqrt(sd) sd=sd/two**ndzero se=se/two**ndzero 902 zmean=sx/zrep+cc zmean=zmean/two**ndzero 910 write(8,*)'no of pops is ',k,' seed is',jseed if(ndenom.gt.0)then write(8,*)' df for var est is',ndenom else write(8,*)' variance assumed known' end if write(8,*)'no of ran dir. is', mocar write(8,*)'value of integral is calculated', irep,' times' write(8,*)'mean value is ',zmean if(zrep.eq.zone)go to 940 write(8,*)'standard error of mean is ',se 940 write(8,*)'--------------------' 9753 format(1x,f3.0,f4.0,i3,f9.6,f9.7,f10.3) 950 continue C WE ARE NOW FINISHED AND NEED TO GO BACK for new instructions deallocate(gd,dc,cmc,dmd,qans,bound) GO TO 21 997 CLOSE(8) CLOSE(9) 998 stop 8000 format(10(f6.4,1x)) end C**************************************** Subroutine pf(n1,n2,z,poff,ier) C (This is a routine to obtain probabilities for the F-distribution) 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 Subroutine to generate direction cosines of a random direction in k C dimensions subroutine racos use comm real sum 10 sum=0. do 20 i=1,k dc(i)=rnor(nz,nw) sum=sum+dc(i)*dc(i) 20 continue sum=sqrt(sum) do i=1,k dc(i)=dc(i)/sum end do 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 20 i=1,i1 z=c1+i x=alog(z) gln=gln+x 20 continue 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 ****************************************************************************** subroutine chol C decomposes pos def matrix vcc into c times c transpose use comm real,dimension(k)::dd integer i,ii,j cmc=0. do 20 i=1,k cmc(i,i)=1. 20 continue do 100 j=1,k dd(j)=vcc(j,j) do 50 ii=1,j-1 dd(j)=dd(j)-cmc(j,ii)*cmc(j,ii)*dd(ii) 50 continue do 80 i=j+1,k cmc(i,j)=vcc(i,j) do 60 ii=1,j-1 cmc(i,j)=cmc(i,j)-cmc(j,ii)*cmc(i,ii)*dd(ii) 60 continue cmc(i,j)=cmc(i,j)/dd(j) 80 continue 100 continue do 120 j=1,k a=sqrt(dd(j)) do 110 i=j,k cmc(i,j)=a*cmc(i,j) 110 continue 120 continue C write(8,*)'cholesky decomposition' C write(8,*)'written from chol' C do i=1,k C write(8,*)(cmc(i,j),j=1,k) C end do return end C***************************************************************** subroutine mvidist(dis1,dis2) C This subroutine calculates the reciprocal distances from the origin C to the boundaries of a convex region with origin wholly within. use comm real dis1,dis2,a integer i,j dis1=0. dis2=0. do 40 i=1,mm a=0. do 30 j=1,k a=a+dc(j)*gd(i,j) 30 continue a=a/dmd(i) if(a.lt.0.) then dis1=min(a,dis1) else dis2=max(a,dis2) endif 40 continue dis1=-dis1 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(8,*) 'input error 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(8,*)'input error 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(8,*)'input error 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(8,*)'A too large, ITMAX too small in GSER' 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(8,*)'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 subroutine crude C This subroutine is used when the convex region does not include the C origin. Instead of direction cosines, a random point is used. C It may also be used as the method of choice use comm real zin,ranpnt,zoneneg integer ifirst zoneneg=-zone zin=zero do ijk=1,mocar if(ndenom.lt.0) then call normpt else call tpoint2(nz,nw,nran,dc,zn,k) end if ifirst=0 80 do iee=1,negs if(ranpnt(iee).gt.zoneneg)go to 100 end do do iee=negs+1,negs+ndzero if(ranpnt(iee).gt.zero)go to 100 end do do iee=negs+ndzero+1,mm if(ranpnt(iee).gt.zone)go to 100 end do C Passed all checks, so "in" zin=zin+zone C By symmetry, cannot have "in" with both dc and -dc when negs.gt.0 C or ndzero.gt.0 if((negs.gt.0).or.(ndzero.gt.0))go to 300 100 if(ifirst.gt.0) go to 300 ifirst=1 dc=-dc go to 80 300 end do points=two*mocar prob=zin/points call date_and_time(date,time) call system_clock(count) ncount2=count seconds=ncount2-ncount1 if(seconds.lt.0)seconds=seconds+count_max seconds=seconds/count_rate write(8,8989)elapse,seconds 8989 format(a28,f10.3) write(8,*)'no of pops is',k,' seed is',jseed if(ndenom.gt.0)then write(8,*)'df for var est is ',ndenom else write(8,*)'variance assumed known' end if write(8,*)'no of random directions is ',mocar write(8,*)'points = ',points write(8,*)'estimate of integral is ',prob sd=sqrt(prob*(zone-prob)/points) write(8,*)'standard error is ',sd write(8,*)'--------------------' return end function ranpnt(iee) C Function is used in crude to determine if point is in region use comm C ranpnt=sum(gd(iee,:)*dc(:)) ranpnt=sum(bound(iee,:)*dc(:)) return end subroutine normpt use comm do i=1,k dc(i)=rnor(nz,nw) end do return end C subroutine tpoint2 to obtain a random t vector for crudel subroutine tpoint2(iz,iw,jran,dc,df,k) real dc(k),df,ch,two integer iz,iw,jran,k two=2. do i=1,k dc(i)=rnor(iz,iw) end do ch=two*gamgen(df/two,iz,iw,jran) ch=sqrt(ch/df) dc=dc/ch return end subroutine binning use comm C Find q, distance from origin of nearest plane q=999999 do i=1,mm if(dmd(i).lt.q)q=dmd(i) end do C Use empirical formula to obtain nq, # of bins for quadrature C NQ=-14.9-39.4*log10(q)+21.8*(log10(q))**2 C Following 4 statements added so that user may elect to set their C own values for nq (say xy)and not use empirical formula. C To do, change 'nqfix=1' to 'nqfix=2' and insert new statement C 'nq=xy'. nqfix=1 if(nqfix.ne.1)go to 100 znq=-14.9-39.4*log10(q)+21.8*(log10(q))**2 if(znq.lt.25)then nq=25 else if(znq.lt.151.3)then C <151.3 corresponds to distance from origin of < .01 nq=znq+10. else mmmvib=abs(mmmvib) go to 910 end if 100 write(8,*)'Binning method used' qq=1/q C Find the probability within the sphere xxx=zm/2 if(ndenom.gt.0) then qsphere = q*q/zm call pf(k,ndenom,qsphere,poff,ier) else qsphere=q*q/two if(qsphere.ge.xxx) then poff=gammp(xxx,qsphere) else poff=zone-gammq(xxx,qsphere) endif end if C q=q/two allocate(wl(nq),pl(nq),prop(nq),ibin(nq)) do 900 njj=1,irep C Set weights and locations call gauleg(zero,qq,pl,wl,nq) C Note that gauleg is used here only in connection with setting C up the bins. C4439 write(8,*)'locations' C4445 write(8,*)(pl(i),i=1,nq) C4446 write(8,*)'weights' C4447 write(8,*)(wl(i),i=1,nq) C Zero bin locations 590 call setbins(nq) do 600 i=1,mocar call racos C write(8,*)'direction cosines are' C write(8,*)(dc(ijk),ijk=1,k) C find reciprocal distance to nearest face call mvidist(dis1,dis2) call fillbin2(dis1,nq) call fillbin2(dis2,nq) 600 continue points=mocar*2 CX zeroprop=zeropts/pts C4501 write(8,*)'bin counts' C4502 write(8,*)(ibin(i),i=1,nq) C4577 write(8,*)'points= ',points,'zeropoints = ',zeropts C Calculate bin proportions call binprop(points,nq) C4509 write(8,*)'bin proportions' C4510 write(8,*)(prop(i),i=1,nq) C call gauleg(zero,qq,pl,wl,nq) C x=fringe(nq) if(ndenom.ge.0) then x=fringe(nq) else x=fringinf(nq) endif CX x=(zone-poff)*zeroprop + x*(zone-zeroprop) qans(njj)=x+poff cc=qans(1) cm=qans(njj)-cc sx=sx+cm ssx=ssx+cm*cm 900 continue deallocate(wl,pl,prop,ibin) 910 return end C This routine zeros the bins for reciprocal distances subroutine setbins(nr) use comm integer i do 20 i=1,nr ibin(i)=0 20 continue return end C**************************************** C This routine does the binning of the reciprocal distances subroutine fillbin2(recdist,nr) use comm integer i real recdist C write(6,*)'recdist = ', recdist do 40 i=nr,1,-1 if(recdist .lt. pl(i)) go to 40 j=nq+1-i ibin(j)=ibin(j)+1 C write(8,*)'in fillbin2,j= ',j,'i= ',i,'pl(i)= ',pl(i) go to 50 40 continue 50 return end C**************************************** C This routine calculates proportions in the bins subroutine binprop(points,nr) C COMMON/GAU/WL(30),PROP(30),PL(30) use comm real points integer i do 10 i=2,nr ibin(i)=ibin(i-1)+ibin(i) 10 continue do 20 i=1,nr prop(i)=1.-float(ibin(i))/points 20 continue return end C**************************************** C This routine calculates the probability of correct statements C outside of region of radius q. function fringe(nr) use comm real x,xx,cd,ce integer j zone=1. two=2. x=0. do 20 i=1,nr j=nr-i+1 ce=zm*alog(zone+zn*pl(i)*pl(i))/two cd=zn*alog(zone+zone/zn/pl(i)/pl(i))/two xx=exp(-alog(pl(i))-ce-cd) tx=xx*wl(i)*prop(j) x=x+tx 20 continue fringe=x*fconst return end C******************************** Subroutine gauleg(x1,x2,x,w,n) C From "Numerical Methods" p 125. Gives constants for Gauss-Legendre. 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) 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 function fringinf(nr) use comm double precision x,xx integer i,j x=0d0 do i=1,nr j=nr-i+1 ce=-gln(zmd2)-(zmd2-zone)*alog(two)-(zm+zone)*alog(pl(i)) &-zone/pl(i)/pl(i)/two xx= exp(ce) x=x+xx*wl(i)*prop(j) end do fringinf=x return end subroutine sortary(n,kk,ra,rb) C Sorts ra in ascending order while making the corresponding arrangement C of the array rb(n,kk) real ra(n),rb(n,kk),rrb(kk) integer n,kk l=n/2+1 ir=n 10 continue if(l.gt.1)then l=l-1 rra=ra(l) rrb(1:kk)=rb(l,1:kk) else rra=ra(ir) rrb(1:kk)=rb(ir,1:kk) ra(ir)=ra(1) rb(ir,1:kk)=rb(1,1:kk) ir=ir-1 if(ir.eq.1)then ra(1)=rra rb(1,1:kk)=rrb(1:kk) return end if end if i=l j=l+l 20 if(j.le.ir)then if(j.lt.ir)then if(ra(j).lt.ra(j+1))j=j+1 end if if(rra.lt.ra(j))then ra(i)=ra(j) rb(i,1:kk)=rb(j,1:kk) i=j j=j+j else j=ir+1 end if go to 20 end if ra(i)=rra rb(i,1:kk)=rrb(1:kk) go to 10 end C gamma generator courtesy of Marsaglia function gamgen(a,z,w,jran) integer z,w,jran real a d=a-.3333333 c=.3333333/sqrt(d) 1 x=rnor(z,w) v=1.+c*x if(v .le. 0.) go to 1 v=v*v*v gamgen=d*v jran=69069*jran U=.5+jran*.2328306e-9 if(U.le.1.-.0331*x**4)return if(log(U)-.5*x*x .le. d*(1.-v+log(v))) return go to 1 return end C Marsaglia's MONTE PYTHON random normal generator 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