c This subroutine implements the Gibbs sampler based on the Polya c urn scheme for the simplest c situation: a Gamma(a,b) mixture of Dirichlet processes whose c alpha measures are exponential(theta) distributions. The total c mass alpha(R) is assumed constant, that is, the mass does c not depend on the value of theta. c c subroutine gibbs1(lo,hi,inlo,inhi,nn,aa0,bb0,totm0, $ nrep,g0,g1,gextra,cc,intim,nc,time,ntim,new, $ surv,sigma,muth,varth, $ doci,tol,guess,beta) c c Declarations for subroutine arguments. c implicit none integer nn,nc,ntim integer nrep,inlo(nn),inhi(nn),intim(nc),g0,g1,gextra double precision lo(nn),hi(nn),aa0,bb0, $ totm0,cc(nc),time(ntim+1), $ surv(ntim),sigma(ntim),muth,varth, $ tol,guess(2,ntim),beta(6,ntim) logical new,doci c c *** INPUT VALUES *** c lo, hi: Vectors containing the interval information. c nn: the number of data observations. c Observation i lies in the interval (lo(i), hi(i)]. c aa0, bb0: the parameters of the Gamma prior for theta. c totm0: the fixed value of alpha(R). (totm0 is short for total mass.) c nrep: the number of simulation repetitions desired. c new: If new is .true., we assume this is a new gibbs run and c the program makes up a starting state. If new is .false., then c the program continues from the last state it reached. c c *** OUTPUT VALUES *** c time: the distinct endpoints of the intervals in the data c plus the additional time points in adtim. These are c ordered from smallest to largest. c ntim: the length of the array time. c surv: the posterior mean of the survival function (at each c of the time points in the array time). c sigma: the posterior standard deviation of the survival function. c mu, var: the posterior mean and variance of theta. c c c save integer nmax,dmx,resumd,mxtim double precision infty,tie parameter(nmax=210,dmx=2*nmax,infty= -99.0,tie= -33.0, $ resumd=40,mxtim=3*nmax) integer ncm integer amass(nmax+1),anxt(nmax+1),aprv(nmax+1),apos(nmax), $ cmass(dmx),cbeg(dmx),cend(dmx),cpos(nmax+1), $ i,j,jrep,jj,k,ii,iimass,ap,cp,numdis,nsumd, $ zbeg,zend,g,isect(nmax+1,2),acnt(nmax+1),cpold,kk, $ bot,top,tp,tbeg logical first,dosrv double precision aval(nmax+1),mcont,theta,u,cum, $ prob(mxtim),tmpsrv(mxtim),am,am2,intp,cump,totp, $ den,den2,sumdis, $ tmp(6),mutot,mu(mxtim) c c Declaring types for functions used. c integer cindex double precision rexpi,pexpi,rgam real ranf c c dosrv=.not.doci if(new)then if (cc(nc) .eq. infty) then ncm=nc-1 else ncm=nc endif den=(nn+totm0)*nrep den2=den*(nn+totm0+1) nsumd=0 mutot=nn+totm0 if (time(1).eq.0.0) then tbeg=2 else tbeg=1 endif c c Create "fake" complete data to be the starting point c for the Gibbs sampler. c All the data values will be distinct (or at least c treated as if they are distinct) and lie c in the appropriate intervals. c Just take value to be the midpoint of the interval c or add one to the lower endpoint if upper endpoint is c infinite. c do 20 k=1,nc cmass(k)=0 20 continue sumdis=0.0d0 numdis=nn do 111 i=1,nn if (hi(i) .eq. infty) then aval(i)=lo(i)+1.0 else aval(i)=(lo(i)+hi(i))/2.0 endif sumdis=sumdis+aval(i) amass(i)=1 acnt(i)=0 cp=cindex(aval(i),inlo(i),inhi(i),cc,nc) apos(i)=i cpos(i)=cp if(cmass(cp) .eq. 0) then cbeg(cp)=i cend(cp)=i else aprv(i)=cend(cp) anxt(cend(cp))=i cend(cp)=i endif cmass(cp)=cmass(cp)+1 111 continue amass(nn+1)=0 acnt(nn+1)=0 zbeg=nn+1 zend=nn+1 endif c c The following quantities are initialized to zero every c time this subroutine is called, even when new=.false.. c if (dosrv) then muth=0.0d0 varth=0.0d0 do 3 k=1,ntim surv(k)=0.0d0 sigma(k)=0.0d0 3 continue else do 333 k=1,ntim do 444 j=1,6 beta(j,k)=0.0d0 444 continue 333 continue endif c c c Begin Gibbs iterations (sweeps). c The iterations 1 to g0-1 are discarded. This is the warm-up. c After that the output statistics are recorded (updated) every c g1 iterations. So g0 is the warm-up and g1 is the spacing. c Every gextra iterations of the "basic" gibbs sampler, we do the c "extra step". A single iteration of the "complete" Gibbs c will then involve gextra iterations of the basic Gibbs sampler c followed by an "extra step". After every g1 iterations of the c complete sampler, we record the output statistics. This means c that output statistics are recorded every g1*gextra iterations c of the basic sampler. c first=.true. g=g0 do 900 jrep=1,nrep do 890 jj=1,g do 880 kk=1,gextra c c Generate theta. c theta=rgam(aa0+numdis,bb0+sumdis) c c Generate the "complete" data values. c do 200 ii=1,nn c c Delete the observation. c ap=apos(ii) cp=cpos(ap) amass(ap)=amass(ap)-1 cmass(cp)=cmass(cp)-1 if(amass(ap) .eq. 0) then sumdis=sumdis-aval(ap) numdis=numdis-1 if(cmass(cp) .ge. 1) then if(cbeg(cp) .eq. ap) then cbeg(cp)=anxt(ap) else if(cend(cp) .eq. ap) then cend(cp)=aprv(ap) else anxt(aprv(ap))=anxt(ap) aprv(anxt(ap))=aprv(ap) endif endif anxt(zend)=ap zend=ap endif c c Generate a new value for this observation. c iimass=0 do 50 k=inlo(ii),inhi(ii) iimass=iimass+cmass(k) 50 continue mcont=totm0*pexpi(theta,lo(ii),hi(ii)) if (ranf() .le. mcont/(mcont+iimass)) then ap=zbeg zbeg=anxt(zbeg) aval(ap)=rexpi(theta,lo(ii),hi(ii)) amass(ap)=1 cp=cindex(aval(ap),inlo(ii),inhi(ii),cc,nc) cpos(ap)=cp if (cmass(cp).eq.0) then cbeg(cp)=ap else anxt(cend(cp))=ap aprv(ap)=cend(cp) endif cend(cp)=ap sumdis=sumdis+aval(ap) numdis=numdis+1 else u=ranf()*iimass cum=0.0 do 60 k=inlo(ii),inhi(ii) cum=cum+cmass(k) if (u .le. cum) then cp=k goto 61 endif 60 continue 61 continue u=ranf()*cmass(cp) ap=cbeg(cp) cum=amass(ap) 69 if(u .gt. cum) then ap=anxt(ap) cum=cum+amass(ap) goto 69 endif amass(ap)=amass(ap)+1 endif apos(ii)=ap cmass(cp)=cmass(cp)+1 200 continue c c Every resumd iterations recompute the value of sumdis. c This guards against accumulation of round-off error. c nsumd=nsumd+1 if (nsumd .eq. resumd) then nsumd=0 sumdis=0.0d0 do 747 ap=1,nn+1 if (amass(ap).ge.1) then sumdis=sumdis+aval(ap) endif 747 continue endif c 880 continue c c Now add the "extra step". Re-randomize the positions c of the atoms, respecting the current group structure. c do 2900 ii=1,nn ap=apos(ii) if(acnt(ap).eq.0)then isect(ap,1)=inlo(ii) isect(ap,2)=inhi(ii) else isect(ap,1)=max(isect(ap,1),inlo(ii)) isect(ap,2)=min(isect(ap,2),inhi(ii)) endif acnt(ap)=acnt(ap)+1 if(acnt(ap).eq.amass(ap))then acnt(ap)=0 sumdis=sumdis-aval(ap) aval(ap)=rexpi(theta,cc(isect(ap,1)),cc(isect(ap,2)+1)) sumdis=sumdis+aval(ap) if(isect(ap,2).gt.isect(ap,1))then cpold=cpos(ap) cp=cindex(aval(ap),isect(ap,1),isect(ap,2),cc,nc) if(cp.ne.cpold)then cpos(ap)=cp cmass(cpold)=cmass(cpold)-amass(ap) if(cmass(cpold) .ge. 1) then if(cbeg(cpold) .eq. ap) then cbeg(cpold)=anxt(ap) else if(cend(cpold) .eq. ap) then cend(cpold)=aprv(ap) else anxt(aprv(ap))=anxt(ap) aprv(anxt(ap))=aprv(ap) endif endif if (cmass(cp).eq.0) then cbeg(cp)=ap else anxt(cend(cp))=ap aprv(ap)=cend(cp) endif cend(cp)=ap cmass(cp)=cmass(cp)+amass(ap) endif endif endif 2900 continue c 890 continue c if (first .eqv. .true.) then g=g1 first=.false. endif if (dosrv) then c c Accumulate survival function info and other output values. c This uses "Fancy" Rao-Blackwellization which involves c conditioning on the group structure. c muth=muth+theta varth=varth+theta*theta c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c prob(1)=exp(-theta*time(1)) tmpsrv(1)=0.0 do 650 i=2,ntim tmpsrv(i)=0.0 prob(i)=exp(-theta*time(i)) prob(i-1)=prob(i-1)-prob(i) 650 continue do 700 ap=1,nn+1 if (amass(ap).ge.1) then bot=intim(isect(ap,1)) top=intim(isect(ap,2)+1) if(bot.eq.(top-1))then tmpsrv(bot)=tmpsrv(bot)+amass(ap) c Note: sigma does not have to be incremented in this case. else totp=pexpi(theta,time(bot),time(top)) cump=0.0 am=amass(ap) am2=am*am do 690 i=top-1,bot+1,-1 intp=prob(i)/totp cump=cump+intp sigma(i)=sigma(i)+am2*cump*(1-cump) tmpsrv(i)=tmpsrv(i)+am*intp 690 continue tmpsrv(bot)=tmpsrv(bot)+am*prob(bot)/totp endif endif 700 continue tmpsrv(ntim)=tmpsrv(ntim)+totm0*prob(ntim) surv(ntim)=surv(ntim)+tmpsrv(ntim) sigma(ntim)=sigma(ntim)+tmpsrv(ntim)+tmpsrv(ntim)**2 do 730 i=ntim-1,1,-1 tmpsrv(i)=tmpsrv(i)+tmpsrv(i+1)+totm0*prob(i) surv(i)=surv(i)+tmpsrv(i) sigma(i)=sigma(i)+tmpsrv(i)+tmpsrv(i)**2 730 continue c else c do 6500 i=1,ntim mu(i)=0.0 6500 continue do 7000 ap=1,nn+1 if (amass(ap).ge.1) then cp=cpos(ap) bot=intim(cp) top=intim(cp+1)-1 if(bot.eq.top)then mu(bot)=mu(bot)+amass(ap) else tp=cindex(aval(ap),bot,top,time,ntim) mu(tp)=mu(tp)+amass(ap) endif endif 7000 continue do 7300 i=ntim-1,1,-1 mu(i)=mu(i)+mu(i+1) 7300 continue do 731 i=1,ntim mu(i)=mu(i)+totm0*exp(-theta*time(i)) 731 continue c c Now compute the beta cdf, etc. needed to compute the c probability intervals. c do 821 i=tbeg,ntim call vbeta(guess(1,i),mu(i),mutot-mu(i),tol,tmp) do 822 j=1,6 beta(j,i)=beta(j,i)+tmp(j) 822 continue 821 continue c c End the case where doci = .TRUE. endif c 900 continue if (dosrv) then c c Now convert totals to means (normalize various quantities). c do 910 k=1,ntim surv(k)=surv(k)/den sigma(k)=sigma(k)/den2 c c Warning: The above line does not actually compute the posterior sigma. c We must convert this into the true sigma (in the Splus program). To c compute sigma in this program, use the following two lines of code: c sigma(k)=sigma(k)/den2-surv(k)**2 c sigma(k)=dsqrt(dmax1(0.0d0,sigma(k))) c 910 continue muth=muth/nrep varth=varth/nrep c Warning: This must be converted to var by subtracting mu*mu as c done below. c var=var/nrep - muth*muth c else c Begin case for doci = .TRUE.. c c Now compute the probability intervals. c Well, maybe not. Just normalize the quantities in beta c and pass back to Splus. c do 9100 j=1,6 do 911 k=tbeg,ntim beta(j,k)=beta(j,k)/nrep 911 continue 9100 continue c c End case for doci = .TRUE.. c endif return end cccccccccccccccccccccccccccccccccc c This subroutine implements the SIS algorithm for the simplest c situation: a Gamma(a,b) mixture of Dirichlet processes whose c alpha measures are exponential(theta) distributions. The total c mass alpha(R) is assumed constant, that is, the mass does c not depend on the value of theta. c c subroutine sis1m(lo,hi,inlo,inhi,nn,aa0,bb0,totm0,aa1,bb1, $ totm1,nrep,cc,intim,nc,domcse,time,ntim, $ surv,sigma,semu,ess,pmax,muth,varth, $ doci,tol,guess,beta) c c Declarations for subroutine arguments. c implicit none integer nn,nc,ntim integer nrep,inlo(nn),inhi(nn),intim(nc) double precision lo(nn),hi(nn),aa0,bb0,aa1,bb1, $ totm0,totm1,cc(nc),time(ntim+1) double precision surv(ntim),sigma(ntim),semu(ntim), $ ess,pmax,muth,varth, $ tol,guess(2,ntim),beta(6,ntim) logical domcse,doci c c *** INPUT VALUES *** c lo, hi: Vectors containing the interval information. c nn: the number of data observations. c Observation i lies in the interval (lo(i), hi(i)]. c Observation i lies in the interval (lo(i), hi(i)]. c cc is a vector containing the distinct values which c appear in lo and hi ordered from smallest to largest. c nc is the length of cc (the number of distinct endpoints). c ncm is the number of distinct endpoints excluding infinity. c inlo(i) and inhi(i) are vectors which give the positions of c lo(i) and hi(i) within cc. c More precisely, inhi(i)+1 is the position of hi(i) in c cc. Then (lo(i),hi(i)] is the union of (cc(k),cc(k+1)] c for k=inlo(i) to inhi(i) c aa0, bb0: the parameters of the Gamma prior for theta. c aa1, bb1: the parameters of the Gamma proposal distribution. c totm0: the fixed value of alpha(R). (totm is short for total mass.) c (this is for the prior distribution). c totm1: the value of alpha(R) used in the proposal distribution. c nrep: the number of simulation repetitions desired. c nc: the "physical" length of the arrays cc and surv. c c *** OUTPUT VALUES *** c This version has been specialized so that the output gives only c surv: the estimated survival function of a "future" observation c computed at the interval endpoints in the data set. c ess: the effective sample size. c pmax: the largest importance sampling probability. c muth, varth: the posterior mean and variance of theta. c c c integer nmax,dmx,mxtim double precision infty,tie parameter(nmax=210,dmx=2*nmax,mxtim=3*nmax, $ infty= -99.0,tie= -33.0) integer amass(nmax),anxt(nmax),cmass(dmx),cm2(dmx),cbeg(dmx), $ cend(dmx),apos,ap,cp,i,jrep,ck,k,ii,iimass,a,tail, $ bot,top,ncm,j,cpos(nmax),tp,cm,tbeg double precision aval(nmax),mcont0,mcont1,pcont0,pcont1, $ theta,u,cum,pei double precision vv,totvv,totv2,den,den2, $ prob(mxtim),tmpsrv(mxtim),intp,cump,totp,tmp,tmp2, $ v2s1(mxtim),v2s2(mxtim), $ btmp(6),mutot,mu(mxtim) logical dosrv c c Declaring types for functions used. c integer cindex double precision rexpi,pexpi,rgam real ranf double precision drat c c dosrv=.not.doci if (cc(nc) .eq. infty) then ncm=nc-1 else ncm=nc endif c c Initialize quantities which will be accumulated c during the SIS iterations. c totvv=0.0d0 totv2=0.0d0 pmax=0.0d0 if (dosrv) then muth=0.0d0 varth=0.0d0 do 3 k=1,ntim surv(k)=0.0d0 sigma(k)=0.0d0 3 continue if(domcse)then do 31 i=1,ntim semu(i)=0.0d0 v2s2(i)=0.0d0 v2s1(i)=0.0d0 31 continue endif else do 333 k=1,ntim do 444 j=1,6 beta(j,k)=0.0d0 444 continue 333 continue mutot=nn+totm0 if (time(1).eq.0.0) then tbeg=2 else tbeg=1 endif c endif c c Begin SIS iterations. c do 900 jrep=1,nrep c c do 10 i=1,nn amass(i)=0 anxt(i)=0 aval(i)=0.0 10 continue do 20 k=1,ncm cmass(k)=0 cm2(k)=0 cbeg(k)=0 cend(k)=0 20 continue c c Generate theta. c theta=rgam(aa1,bb1) vv=drat(theta,aa0,bb0,aa1,bb1) c c Generate first data value. c c (Note: The first data value is always generated c from the continuous part of the urn, so c there is no need to include an additional c factor involving totm1 versus totm0.) c aval(1)=rexpi(theta,lo(1),hi(1)) vv=vv*pexpi(theta,lo(1),hi(1)) amass(1)=1 ck=cindex(aval(1),inlo(1),inhi(1),cc,nc) c The next statement is needed only for ci's. cpos(1)=ck cmass(ck)=1 cm2(ck)=1 cbeg(ck)=1 cend(ck)=1 c c Generate other data values. c apos=1 do 200 ii=2,nn iimass=0 do 50 k=inlo(ii),inhi(ii) iimass=iimass+cmass(k) 50 continue pei=pexpi(theta,lo(ii),hi(ii)) mcont0=totm0*pei mcont1=totm1*pei pcont0=mcont0/(mcont0+iimass) pcont1=mcont1/(mcont1+iimass) vv=vv*(mcont0+iimass)/(totm0+ii-1) if (ranf() .le. pcont1) then apos=apos+1 vv=vv*(pcont0/pcont1) aval(apos)=rexpi(theta,lo(ii),hi(ii)) amass(apos)=1 ck=cindex(aval(apos),inlo(ii),inhi(ii),cc,nc) c The next statement is needed only for ci's. cpos(apos)=ck cm2(ck)=cm2(ck)+1 if (cmass(ck).eq.0) then cbeg(ck)=apos else anxt(cend(ck))=apos endif cend(ck)=apos else vv=vv*(1-pcont0)/(1-pcont1) u=ranf()*iimass cum=0.0 do 60 k=inlo(ii),inhi(ii) cum=cum+cmass(k) if (u .le. cum) then ck=k goto 61 endif 60 continue 61 continue u=ranf()*cmass(ck) a=cbeg(ck) cum=amass(a) 69 if(u .gt. cum) then a=anxt(a) cum=cum+amass(a) goto 69 endif amass(a)=amass(a)+1 c Accumulate sum of squared atom masses in each interval. cm2(ck)=cm2(ck)+2*amass(a)-1 endif cmass(ck)=cmass(ck)+1 200 continue c c Accumulate survival function info and other output values. c totvv=totvv+vv totv2=totv2+vv*vv pmax=dmax1(pmax,vv) c if (dosrv) then muth=muth+vv*theta varth=varth+vv*theta*theta tail=0 c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c prob(1)=exp(-theta*time(1)) tmpsrv(1)=0.0 do 650 i=2,ntim tmpsrv(i)=0.0 prob(i)=exp(-theta*time(i)) prob(i-1)=prob(i-1)-prob(i) 650 continue do 652 ck=1,ncm if(cmass(ck).ge.1)then cm=cmass(ck) bot=intim(ck) top=intim(ck+1) if(bot.eq.(top-1))then tmpsrv(bot)=tmpsrv(bot)+cm c Note: sigma does not have to be incremented in this case. else totp=pexpi(theta,time(bot),time(top)) cump=0.0 do 690 i=top-1,bot+1,-1 intp=prob(i)/totp cump=cump+intp sigma(i)=sigma(i)+vv*cm2(ck)*cump*(1-cump) tmpsrv(i)=tmpsrv(i)+cm*intp 690 continue tmpsrv(bot)=tmpsrv(bot)+cm*prob(bot)/totp endif endif 652 continue tmpsrv(ntim)=tmpsrv(ntim)+totm0*prob(ntim) surv(ntim)=surv(ntim)+vv*tmpsrv(ntim) sigma(ntim)=sigma(ntim)+ $ vv*(tmpsrv(ntim)+tmpsrv(ntim)**2) do 730 i=ntim-1,1,-1 tmpsrv(i)=tmpsrv(i)+tmpsrv(i+1)+totm0*prob(i) surv(i)=surv(i)+vv*tmpsrv(i) sigma(i)=sigma(i)+vv*(tmpsrv(i)+tmpsrv(i)**2) 730 continue if(domcse) then do 735 i=1,ntim tmp=vv*tmpsrv(i) v2s1(i)=v2s1(i)+vv*tmp v2s2(i)=v2s2(i)+tmp*tmp 735 continue endif c else c do 6500 i=1,ntim mu(i)=0.0 6500 continue do 7000 ap=1,apos cp=cpos(ap) bot=intim(cp) top=intim(cp+1)-1 if(bot.eq.top)then mu(bot)=mu(bot)+amass(ap) else tp=cindex(aval(ap),bot,top,time,ntim) mu(tp)=mu(tp)+amass(ap) endif 7000 continue do 7300 i=ntim-1,1,-1 mu(i)=mu(i)+mu(i+1) 7300 continue do 731 i=1,ntim mu(i)=mu(i)+totm0*exp(-theta*time(i)) 731 continue c c Now compute the beta cdf, etc. needed to compute the c probability intervals. c do 821 i=tbeg,ntim call vbeta(guess(1,i),mu(i),mutot-mu(i),tol,btmp) do 822 j=1,6 beta(j,i)=beta(j,i)+vv*btmp(j) 822 continue 821 continue c c End the case where doci = .TRUE. endif c c End of SIS iteration. Now start another. c 900 continue c c Now divide through by totvv and normalize various quantities. c pmax=pmax/totvv ess=totvv**2/totv2 if (dosrv) then muth=muth/totvv varth=varth/totvv - muth*muth den=nn+totm0 den2=den*(den+1) do 950 k=1,ntim surv(k)=surv(k)/(totvv*den) sigma(k)=sigma(k)/(totvv*den2) - surv(k)**2 sigma(k)=dsqrt(dmax1(0.0d0,sigma(k))) 950 continue if(domcse)then tmp=den*den tmp2=totvv*totvv do 955 i=1,ntim v2s1(i)=v2s1(i)/den v2s2(i)=v2s2(i)/tmp semu(i)=v2s2(i) + (surv(i)**2)*totv2 $ - 2.0*surv(i)*v2s1(i) semu(i)=dsqrt(dmax1(0.0d0,semu(i)/tmp2)) 955 continue endif else c do 9100 j=1,6 do 911 k=tbeg,ntim beta(j,k)=beta(j,k)/totvv 911 continue 9100 continue c c End case for doci = .TRUE.. c endif c return end cccccccccccccccccccccccccccccccccc c This subroutine implements the SIS algorithm for the simplest c situation: a Gamma(a,b) mixture of Dirichlet processes whose c alpha measures are exponential(theta) distributions. The total c mass alpha(R) is assumed constant, that is, the mass does c not depend on the value of theta. c c subroutine sis1(lo,hi,inlo,inhi,nn,aa0,bb0,totm0, $ nrep,cc,intim,nc,domcse,time,ntim, $ surv,sigma,semu,ess,pmax,muth,varth, $ doci,tol,guess,beta) c c Declarations for subroutine arguments. c implicit none integer nn,nc,ntim integer nrep,inlo(nn),inhi(nn),intim(nc) double precision lo(nn),hi(nn),aa0,bb0, $ totm0,cc(nc),time(ntim+1) double precision surv(ntim),sigma(ntim),semu(ntim), $ ess,pmax,muth,varth, $ tol,guess(2,ntim),beta(6,ntim) logical domcse,doci c c *** INPUT VALUES *** c lo, hi: Vectors containing the interval information. c nn: the number of data observations. c Observation i lies in the interval (lo(i), hi(i)]. c Observation i lies in the interval (lo(i), hi(i)]. c cc is a vector containing the distinct values which c appear in lo and hi ordered from smallest to largest. c nc is the length of cc (the number of distinct endpoints). c ncm is the number of distinct endpoints excluding infinity. c inlo(i) and inhi(i) are vectors which give the positions of c lo(i) and hi(i) within cc. c More precisely, inhi(i)+1 is the position of hi(i) in c cc. Then (lo(i),hi(i)] is the union of (cc(k),cc(k+1)] c for k=inlo(i) to inhi(i) c aa0, bb0: the parameters of the Gamma prior for theta. c totm0: the fixed value of alpha(R). (totm is short for total mass.) c (this is for the prior distribution). c nrep: the number of simulation repetitions desired. c nc: the "physical" length of the arrays cc and surv. c c *** OUTPUT VALUES *** c This version has been specialized so that the output gives only c surv: the estimated survival function of a "future" observation c computed at the interval endpoints in the data set. c ess: the effective sample size. c pmax: the largest importance sampling probability. c muth, varth: the posterior mean and variance of theta. c c c integer nmax,dmx,mxtim double precision infty,tie parameter(nmax=210,dmx=2*nmax,mxtim=3*nmax, $ infty= -99.0,tie= -33.0) integer amass(nmax),anxt(nmax),cmass(dmx),cm2(dmx),cbeg(dmx), $ cend(dmx),apos,ap,cp,i,jrep,ck,k,ii,iimass,a,tail, $ bot,top,ncm,j,cpos(nmax),tp,cm,tbeg double precision aval(nmax),mcont0,pcont0, $ theta,u,cum,pei double precision vv,totvv,totv2,den,den2, $ prob(mxtim),tmpsrv(mxtim),intp,cump,totp,tmp,tmp2, $ v2s1(mxtim),v2s2(mxtim), $ btmp(6),mutot,mu(mxtim) logical dosrv c c Declaring types for functions used. c integer cindex double precision rexpi,pexpi,rgam real ranf c c dosrv=.not.doci if (cc(nc) .eq. infty) then ncm=nc-1 else ncm=nc endif c c Initialize quantities which will be accumulated c during the SIS iterations. c totvv=0.0d0 totv2=0.0d0 pmax=0.0d0 if (dosrv) then muth=0.0d0 varth=0.0d0 do 3 k=1,ntim surv(k)=0.0d0 sigma(k)=0.0d0 3 continue if(domcse)then do 31 i=1,ntim semu(i)=0.0d0 v2s2(i)=0.0d0 v2s1(i)=0.0d0 31 continue endif else do 333 k=1,ntim do 444 j=1,6 beta(j,k)=0.0d0 444 continue 333 continue mutot=nn+totm0 if (time(1).eq.0.0) then tbeg=2 else tbeg=1 endif c endif c c Begin SIS iterations. c do 900 jrep=1,nrep c c do 10 i=1,nn amass(i)=0 anxt(i)=0 aval(i)=0.0 10 continue do 20 k=1,ncm cmass(k)=0 cm2(k)=0 cbeg(k)=0 cend(k)=0 20 continue c c Generate theta. c theta=rgam(aa0,bb0) vv=1.0d0 c c Generate first data value. c c aval(1)=rexpi(theta,lo(1),hi(1)) vv=vv*pexpi(theta,lo(1),hi(1)) amass(1)=1 ck=cindex(aval(1),inlo(1),inhi(1),cc,nc) c The next statement is needed only for ci's. cpos(1)=ck cmass(ck)=1 cm2(ck)=1 cbeg(ck)=1 cend(ck)=1 c c Generate other data values. c apos=1 do 200 ii=2,nn iimass=0 do 50 k=inlo(ii),inhi(ii) iimass=iimass+cmass(k) 50 continue pei=pexpi(theta,lo(ii),hi(ii)) mcont0=totm0*pei pcont0=mcont0/(mcont0+iimass) vv=vv*(mcont0+iimass)/(totm0+ii-1) if (ranf() .le. pcont0) then apos=apos+1 aval(apos)=rexpi(theta,lo(ii),hi(ii)) amass(apos)=1 ck=cindex(aval(apos),inlo(ii),inhi(ii),cc,nc) c The next statement is needed only for ci's. cpos(apos)=ck cm2(ck)=cm2(ck)+1 if (cmass(ck).eq.0) then cbeg(ck)=apos else anxt(cend(ck))=apos endif cend(ck)=apos else u=ranf()*iimass cum=0.0 do 60 k=inlo(ii),inhi(ii) cum=cum+cmass(k) if (u .le. cum) then ck=k goto 61 endif 60 continue 61 continue u=ranf()*cmass(ck) a=cbeg(ck) cum=amass(a) 69 if(u .gt. cum) then a=anxt(a) cum=cum+amass(a) goto 69 endif amass(a)=amass(a)+1 c Accumulate sum of squared atom masses in each interval. cm2(ck)=cm2(ck)+2*amass(a)-1 endif cmass(ck)=cmass(ck)+1 200 continue c c Accumulate survival function info and other output values. c totvv=totvv+vv totv2=totv2+vv*vv pmax=dmax1(pmax,vv) c if (dosrv) then muth=muth+vv*theta varth=varth+vv*theta*theta tail=0 c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c prob(1)=exp(-theta*time(1)) tmpsrv(1)=0.0 do 650 i=2,ntim tmpsrv(i)=0.0 prob(i)=exp(-theta*time(i)) prob(i-1)=prob(i-1)-prob(i) 650 continue do 652 ck=1,ncm if(cmass(ck).ge.1)then cm=cmass(ck) bot=intim(ck) top=intim(ck+1) if(bot.eq.(top-1))then tmpsrv(bot)=tmpsrv(bot)+cm c Note: sigma does not have to be incremented in this case. else totp=pexpi(theta,time(bot),time(top)) cump=0.0 do 690 i=top-1,bot+1,-1 intp=prob(i)/totp cump=cump+intp sigma(i)=sigma(i)+vv*cm2(ck)*cump*(1-cump) tmpsrv(i)=tmpsrv(i)+cm*intp 690 continue tmpsrv(bot)=tmpsrv(bot)+cm*prob(bot)/totp endif endif 652 continue tmpsrv(ntim)=tmpsrv(ntim)+totm0*prob(ntim) surv(ntim)=surv(ntim)+vv*tmpsrv(ntim) sigma(ntim)=sigma(ntim)+ $ vv*(tmpsrv(ntim)+tmpsrv(ntim)**2) do 730 i=ntim-1,1,-1 tmpsrv(i)=tmpsrv(i)+tmpsrv(i+1)+totm0*prob(i) surv(i)=surv(i)+vv*tmpsrv(i) sigma(i)=sigma(i)+vv*(tmpsrv(i)+tmpsrv(i)**2) 730 continue if(domcse) then do 735 i=1,ntim tmp=vv*tmpsrv(i) v2s1(i)=v2s1(i)+vv*tmp v2s2(i)=v2s2(i)+tmp*tmp 735 continue endif c else c do 6500 i=1,ntim mu(i)=0.0 6500 continue do 7000 ap=1,apos cp=cpos(ap) bot=intim(cp) top=intim(cp+1)-1 if(bot.eq.top)then mu(bot)=mu(bot)+amass(ap) else tp=cindex(aval(ap),bot,top,time,ntim) mu(tp)=mu(tp)+amass(ap) endif 7000 continue do 7300 i=ntim-1,1,-1 mu(i)=mu(i)+mu(i+1) 7300 continue do 731 i=1,ntim mu(i)=mu(i)+totm0*exp(-theta*time(i)) 731 continue c c Now compute the beta cdf, etc. needed to compute the c probability intervals. c do 821 i=tbeg,ntim call vbeta(guess(1,i),mu(i),mutot-mu(i),tol,btmp) do 822 j=1,6 beta(j,i)=beta(j,i)+vv*btmp(j) 822 continue 821 continue c c End the case where doci = .TRUE. endif c c End of SIS iteration. Now start another. c 900 continue c c Now divide through by totvv and normalize various quantities. c pmax=pmax/totvv ess=totvv**2/totv2 if (dosrv) then muth=muth/totvv varth=varth/totvv - muth*muth den=nn+totm0 den2=den*(den+1) do 950 k=1,ntim surv(k)=surv(k)/(totvv*den) sigma(k)=sigma(k)/(totvv*den2) - surv(k)**2 sigma(k)=dsqrt(dmax1(0.0d0,sigma(k))) 950 continue if(domcse)then tmp=den*den tmp2=totvv*totvv do 955 i=1,ntim v2s1(i)=v2s1(i)/den v2s2(i)=v2s2(i)/tmp semu(i)=v2s2(i) + (surv(i)**2)*totv2 $ - 2.0*surv(i)*v2s1(i) semu(i)=dsqrt(dmax1(0.0d0,semu(i)/tmp2)) 955 continue endif else c do 9100 j=1,6 do 911 k=tbeg,ntim beta(j,k)=beta(j,k)/totvv 911 continue 9100 continue c c End case for doci = .TRUE.. c endif c return end cccccccccccccccccccccccccccccccccc integer function cindex(x,i1,i2,cc,nc) implicit none integer i1,i2,nc,k double precision x,cc(nc) cindex=i2 do 10 k=i1+1,i2 if (x .le. cc(k)) then cindex=k-1 goto 11 endif 10 continue 11 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function rexpi(t,a,b) implicit none c This function generates a random exponential with c parameter theta conditional on it lying in the interval (a,b). double precision t,a,b double precision infty,tie,u real ranf parameter(infty= -99.0,tie= -33.0) c if (b .eq. infty) then rexpi = a-log(ranf())/t else if (b .gt. a) then u = ranf() rexpi = a-log(u+(1-u)*exp(-t*(b-a)))/t else if (b .eq. tie) then rexpi = a else c write(*,*) "rexpi: bad values for a, b." rexpi=infty endif return end double precision function pexpi(t,a,b) implicit none double precision t,a,b double precision infty,tie parameter(infty= -99.0,tie= -33.0) c if (b .eq. infty) then pexpi = exp(-t*a) else if (b .gt. a) then pexpi = exp(-t*a) - exp(-t*b) else if (b .eq. tie) then pexpi = 0.0 else c write(*,*) "pexpi: bad values for a, b." pexpi=infty endif return end double precision function rgam(a,b) implicit none double precision a,b real sgamma rgam = sgamma(real(a))/b return end double precision function drat(theta,a0,b0,a1,b1) c This computes the log of the ratio of two gamma densities. c Constants not involving theta are dropped. implicit none double precision theta,a0,b0,a1,b1 drat=dexp((a0-a1)*dlog(theta)-theta*(b0-b1)) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This subroutine handles only right-censoring. c It simulates directly from the posterior distribution. c The prior is: a Gamma(a,b) mixture of Dirichlet processes whose c alpha measures are exponential(theta) distributions. The total c mass alpha(R) is assumed constant, that is, the mass does c not depend on the value of theta. c c subroutine ritcen(z,delt,intim,nn,time,ntim, $ aa,bb,totm,nrep,domcse,nthint, $ surv,sigma,semu,sesig,muth,varth,eff, $ doci,tol,guess,beta) c c Declarations for subroutine arguments. c implicit none integer nn,nrep,ntim,intim(nn),delt(nn),nthint double precision z(nn),aa,bb,totm,time(ntim+1), $ surv(ntim),sigma(ntim),semu(ntim),sesig(ntim), $ muth,varth,eff, $ tol,guess(2,ntim),beta(6,ntim) logical domcse,doci c c *** INPUT VALUES *** c z, delt: Data times and indicator of uncensored observations. c nn: the number of data observations. c totm: the fixed value of alpha(R). (totm is short for total mass.) c (this is for the prior distribution). c nrep: the number of simulation repetitions desired. c *** OUTPUT VALUES *** c surv: the estimated survival function of a "future" observation c computed at the various time points. c muth, varth: the posterior mean and variance of theta. c c c save integer nmax,mxtim double precision big parameter(nmax=210,mxtim=3*nmax,big= 2.0**100) integer amass(nmax),apos,i,jrep,ii,a,bot,icen, $ amin(nmax),apin,ncen,atim(nmax),anum(nmax), $ am,am2,cen(nmax),ntry,tp,tbeg,j,k double precision mcont,pcont,theta,prev,x double precision den,den2,den3,den4,den5,tprob(mxtim), $ prob(mxtim),tmpsrv(mxtim),tmpsg2(mxtim), $ intp,cump,totp,cov(mxtim), $ tmp(6),mutot,mu(mxtim) logical dosrv c c Declaring types for functions used. c real rpthet,ranf integer cindex c dosrv=.not.doci if (dosrv) then c c Set up the random number generator for generating from the c posterior distribution of theta. c call rptset(z,delt,nn,aa,bb,totm,nthint) c c Initialize quantities which will be accumulated c during the SIS iterations. c muth=0.0d0 varth=0.0d0 eff=0.0d0 do 30 i=1,ntim surv(i)=0.0d0 sigma(i)=0.0d0 30 continue if(domcse)then do 31 i=1,ntim semu(i)=0.0d0 sesig(i)=0.0d0 cov(i)=0.0d0 31 continue endif c c Set up atoms at the uncensored observations. c Create list of censored observations. c For this code to work properly, observations c must be arranged so that the times are from largest c to smallest, with censored observations placed before c uncensored observations having the same time. (This c corresponds to having censored observations dying c at some point strictly after the censoring time. c apin=0 prev=big ncen=0 do 50 i=1,nn if(delt(i).eq.1)then if(z(i).lt.prev)then prev=z(i) apin=apin+1 amin(apin)=1 c xxx Warning: intim is defined differently than in other programs. atim(apin)=intim(i) else amin(apin)=amin(apin)+1 endif anum(i)=apin else ncen=ncen+1 cen(ncen)=i endif 50 continue else mutot=nn+totm if (time(1).eq.0.0) then tbeg=2 else tbeg=1 endif do 333 k=1,ntim do 444 j=1,6 beta(j,k)=0.0d0 444 continue 333 continue endif c c Begin SIS iterations. c do 900 jrep=1,nrep c theta=rpthet(ntry) c do 70 i=1,apin amass(i)=amin(i) 70 continue c apos=apin do 300 icen=1,ncen ii=cen(icen) mcont=totm*exp(-theta*z(ii)) pcont=mcont/(mcont+ii-1) if(ranf().le.pcont)then apos=apos+1 amass(apos)=1 c xxx intim defined differently. atim(apos)=intim(ii) anum(ii)=apos else a=anum(1+ifix(ranf()*(ii-1))) anum(ii)=a amass(a)=amass(a)+1 endif 300 continue if (dosrv) then c c Accumulate survival function info and other output values. c This uses "Fancy" Rao-Blackwellization which involves c conditioning on the group structure. c eff=eff+ntry muth=muth+theta varth=varth+theta*theta c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c do 401 i=1,ntim tmpsrv(i)=0.0 tmpsg2(i)=0.0 401 continue tprob(1)=exp(-theta*time(1)) do 650 i=2,ntim tprob(i)=exp(-theta*time(i)) prob(i-1)=tprob(i-1)-tprob(i) 650 continue prob(ntim)=tprob(ntim) do 675 a=1,apin bot=atim(a) tmpsrv(bot)=tmpsrv(bot)+amass(a) 675 continue do 700 a=apin+1,apos bot=atim(a) if(bot.eq.ntim)then tmpsrv(bot)=tmpsrv(bot)+amass(a) c Note: tmpsg2 does not have to be incremented in this case. else totp=tprob(bot) cump=0.0 am=amass(a) am2=am*am do 690 i=ntim,bot+1,-1 intp=prob(i)/totp cump=cump+intp tmpsg2(i)=tmpsg2(i)+am2*cump*(1-cump) tmpsrv(i)=tmpsrv(i)+am*intp 690 continue tmpsrv(bot)=tmpsrv(bot)+am*prob(bot)/totp endif 700 continue tmpsrv(ntim)=tmpsrv(ntim)+totm*prob(ntim) tmpsg2(ntim)=tmpsg2(ntim)+tmpsrv(ntim)+tmpsrv(ntim)**2 do 730 i=ntim-1,1,-1 tmpsrv(i)=tmpsrv(i)+tmpsrv(i+1)+totm*prob(i) tmpsg2(i)=tmpsg2(i)+tmpsrv(i)+tmpsrv(i)**2 730 continue do 731 i=1,ntim surv(i)=surv(i)+tmpsrv(i) sigma(i)=sigma(i)+tmpsg2(i) 731 continue if(domcse)then do 732 i=1,ntim semu(i)=semu(i)+tmpsrv(i)**2 sesig(i)=sesig(i)+tmpsg2(i)**2 cov(i)=cov(i)+tmpsrv(i)*tmpsg2(i) 732 continue endif else c c begin doci = .TRUE. c do 6500 i=1,ntim mu(i)=0.0 6500 continue do 6750 a=1,apin bot=atim(a) mu(bot)=mu(bot)+amass(a) 6750 continue do 7000 a=apin+1,apos bot=atim(a) if(bot.eq.ntim)then mu(bot)=mu(bot)+amass(a) else x=time(bot)-log(ranf())/theta tp=cindex(x,bot,ntim,time,ntim) mu(tp)=mu(tp)+amass(a) endif 7000 continue do 7300 i=ntim-1,1,-1 mu(i)=mu(i)+mu(i+1) 7300 continue do 7310 i=1,ntim mu(i)=mu(i)+totm*exp(-theta*time(i)) 7310 continue c c Now compute the beta cdf, etc. needed to compute the c probability intervals. c do 821 i=tbeg,ntim call vbeta(guess(1,i),mu(i),mutot-mu(i),tol,tmp) do 822 j=1,6 beta(j,i)=beta(j,i)+tmp(j) 822 continue 821 continue c c End the case where doci = .TRUE. endif c 900 continue c c Now convert totals to means (normalize various quantities), c variances, and covariance (if domcse=.true.). c if (dosrv) then den=(nn+totm)*nrep den2=(nn+totm)*(nn+totm+1.0d0)*nrep do 910 i=1,ntim surv(i)=surv(i)/den sigma(i)=sigma(i)/den2 910 continue if (domcse) then den3=(den**2)/nrep den4=(den2**2)/nrep den5=den*den2/nrep do 911 i=1,ntim semu(i)=semu(i)/den3 - surv(i)**2 sesig(i)=sesig(i)/den4 - sigma(i)**2 cov(i)=cov(i)/den5 - surv(i)*sigma(i) 911 continue endif c c do 912 i=1,ntim sigma(i)=sigma(i)-surv(i)**2 912 continue if (domcse) then do 913 i=1,ntim sesig(i)=dsqrt(dmax1(0.0d0, $ (sesig(i)-4.0d0*surv(i)*cov(i)+ $ 4.0d0*(surv(i)**2)*semu(i))/(4.0d0*nrep*sigma(i)))) semu(i)=dsqrt(dmax1(0.0d0,semu(i)/nrep)) 913 continue endif do 914 i=1,ntim sigma(i)=dsqrt(dmax1(0.0d0,sigma(i))) 914 continue muth=muth/nrep varth=varth/nrep - muth*muth eff=eff/nrep else c c Begin doci = .TRUE. c do 9100 j=1,6 do 9110 k=tbeg,ntim beta(j,k)=beta(j,k)/nrep 9110 continue 9100 continue c c End doci = .TRUE. c endif return end real function rpthet(n) c *************************** c Sample from the posterior for theta. c *************************** c The argument n records the number of X values we generate c until finally accepting one. This is used to study the c performance of the algorithm. c implicit none integer n,i,ii real ranf,x,llfn,u,pkeep,ginv c c Variables in common. c integer nth,mxi parameter(mxi=50) real cump(mxi),aa,bb(mxi),th(mxi),cc(mxi),dd(mxi), $ ss(mxi),ll(mxi) common /glob1/nth,cump,aa,bb,th,cc,dd,ss,ll save /glob1/ c n=0 10 continue n=n+1 u=ranf() do 30 i=1,nth if(u.le.cump(i))then c .. Generate x from gamma(aa,bb(i)) restricted to c .. the interval (th(i),th(i+1)). x=ginv(aa,cc(i)+dd(i)*ranf())/bb(i) ii=i goto 31 endif 30 continue 31 continue c .. accept or reject pkeep=exp(llfn(x)-(ll(ii)+ss(ii)*(x-th(ii)))) if (ranf().gt.pkeep) goto 10 rpthet=x return end c c real function ginv(a,u) implicit none integer ierr real u,a call gaminv(a,ginv,-1.0,u,1.0-u,ierr) return end c c c c real function gcdf(a,x) implicit none real x,a,q call gratio(a,x,gcdf,q,0) return end c c real function llfn(x) c c This computes the log of the function h which is the c product of terms arising from the censored observations. c implicit none real x integer i c c Variables in common. c integer ncen,mxcen parameter(mxcen=100) c .. The array kk will contain integers, but is declared real. real totm,zz(mxcen),kk(mxcen) common /glob2/ncen,totm,zz,kk save /glob2/ c llfn=0.0 do 10 i=1,ncen llfn=llfn+log(kk(i)+totm*exp(-x*zz(i))) 10 continue return end subroutine rptset(z,delt,nn,aa0,bb0,totm0,nth0) c c z is the array of data times. c delt gives indicators (1 or 0) for the uncensored observations. c Data must be ordered with times z(i) from largest to smallest, c with censored observations placed before uncensored observations c having the same time. c nn is the number of observations. c aa0, bb0, totm0 are the parameters of the prior. c nth0 is the number of theta intervals to be used in the c acceptance/rejection generation scheme. c implicit none integer nn,nth0 double precision z(nn),aa0,bb0,totm0 integer delt(nn),i real big,prev,bb1,lim,sum,llfn,targ,gcdf,mult,tinc c c Variables in common. c integer nth,mxi parameter(mxi=50,big=2.0**30) real cump(mxi),aa,bb(mxi),th(mxi),cc(mxi),dd(mxi), $ ss(mxi),ll(mxi) common /glob1/nth,cump,aa,bb,th,cc,dd,ss,ll save /glob1/ c integer ncen,mxcen parameter(mxcen=100) c .. The array kk will contain integers, but is declared real. real totm,zz(mxcen),kk(mxcen) common /glob2/ncen,totm,zz,kk save /glob2/ c c Give values to variables in glob2, and to aa, bb. c Note, if the largest time is censored, then ncen will c be one less than the number of censored observations. c totm=totm0 aa=aa0 bb1=bb0 ncen=0 prev=big do 100 i=1,nn if(delt(i).eq.1)then if(z(i).lt.prev)then prev=z(i) aa=aa+1 bb1=bb1+z(i) endif else if(i.eq.1)then c .. absorb largest censored observation c .. into gamma dist. bb1=bb1+z(i) else ncen=ncen+1 zz(ncen)=z(i) kk(ncen)=i-1 endif 100 continue c c Give values to remaining variables in glob1. c nth=nth0 c c Determine the array th(i). c lim=0.0 do 150 i=1,ncen lim=lim+log(kk(i)) 150 continue c th(1)=0.0 ll(1)=llfn(0.0) tinc=(ll(1)-lim)/nth do 200 i=2,nth targ=ll(1)-(i-1)*tinc c .. th(i) is the solution to llfn(x)=targ obtained by c .. newton-raphson starting from th(i-1). call llsolv(targ,th(i),th(i-1)) ll(i)=llfn(th(i)) ss(i-1)=(ll(i)-ll(i-1))/(th(i)-th(i-1)) bb(i-1)=bb1-ss(i-1) 200 continue ss(nth)=0.0 bb(nth)=bb1 c sum=0.0 do 250 i=1,nth cc(i)=gcdf(aa,bb(i)*th(i)) if(i.lt.nth)then dd(i)=gcdf(aa,bb(i)*th(i+1))-cc(i) else dd(i)=1.0-cc(i) endif mult=exp(aa*log(bb1/bb(i))+ll(i)-ss(i)*th(i)) sum=sum+mult*dd(i) cump(i)=sum 250 continue do 260 i=1,nth-1 cump(i)=cump(i)/sum 260 continue cump(nth)=1.0 c return end subroutine llsolv(targ,xans,xinit) implicit none real targ,xans,xinit real x,val,llf,lld,y integer i c c Variables in common. c integer ncen,mxcen parameter(mxcen=100) c .. The array kk will contain integers, but is declared real. real totm,zz(mxcen),kk(mxcen) common /glob2/ncen,totm,zz,kk save /glob2/ c x=xinit 5 continue llf=0.0 lld=0.0 do 10 i=1,ncen y=totm*exp(-x*zz(i)) val=kk(i)+y llf=llf+log(val) lld=lld-zz(i)*y/val 10 continue x=x-(llf-targ)/lld if(abs(llf-targ).gt.0.001) goto 5 xans=x return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This subroutine implements the SIS2 algorithm for the simplest c situation: a Gamma(a,b) mixture of Dirichlet processes whose c alpha measures are exponential(theta) distributions. The total c mass alpha(R) is assumed constant, that is, the mass does c not depend on the value of theta. c c subroutine sis2(lo,hi,inlo,inhi,nn,aa0,bb0,totm0, $ nrep,cc,intim,nc,domcse,time,ntim, $ surv,sigma,semu,ess,pmax,muth,varth, $ doci,tol,guess,beta) c c Declarations for subroutine arguments. c implicit none integer nn,nc,ntim integer nrep,inlo(nn),inhi(nn),intim(nc) double precision lo(nn),hi(nn),aa0,bb0, $ totm0,cc(nc),time(ntim+1) double precision surv(ntim),sigma(ntim),semu(ntim), $ ess,pmax,muth,varth, $ tol,guess(2,ntim),beta(6,ntim) logical domcse,doci c c *** INPUT VALUES *** c lo, hi: Vectors containing the interval information. c nn: the number of data observations. c Observation i lies in the interval (lo(i), hi(i)]. c Observation i lies in the interval (lo(i), hi(i)]. c cc is a vector containing the distinct values which c appear in lo and hi ordered from smallest to largest. c nc is the length of cc (the number of distinct endpoints). c ncm is the number of distinct endpoints excluding infinity. c inlo(i) and inhi(i) are vectors which give the positions of c lo(i) and hi(i) within cc. c More precisely, inhi(i)+1 is the position of hi(i) in c cc. Then (lo(i),hi(i)] is the union of (cc(k),cc(k+1)] c for k=inlo(i) to inhi(i) c aa0, bb0: the parameters of the Gamma prior for theta. c totm0: the fixed value of alpha(R). (totm is short for total mass.) c (this is for the prior distribution). c nrep: the number of simulation repetitions desired. c c *** OUTPUT VALUES *** c This version has been specialized so that the output gives only c surv: the estimated survival function of a "future" observation c computed at the interval endpoints in the data set. c ess: the effective sample size. c pmax: the largest importance sampling probability. c mu, var: the posterior mean and variance of theta. c c c integer nmax,dmx,mxtim double precision infty,tie parameter(nmax=210,dmx=2*nmax,mxtim=3*nmax, $ infty= -99.0,tie= -33.0) integer amass(nmax),anxt(nmax),cmass(dmx),cm2(dmx),cbeg(dmx), $ cend(dmx),apos,i,jrep,ck,k,ii,iimass,a,tail,cm, $ bot,top,ncm,j,cpos(nmax),tp,tbeg,cp,ap double precision aval(nmax),mcont,theta,u,cum,aaii,bbii,x double precision vv,totvv,totv2,den,den2, $ prob(mxtim),tmpsrv(mxtim),intp,cump,totp,tmp,tmp2, $ v2s1(mxtim),v2s2(mxtim), $ btmp(6),mutot,mu(mxtim) logical dosrv c c Declaring types for functions used. c integer cindex double precision pexpi,rgam,rpari,ppari real ranf c c dosrv=.not.doci if (cc(nc) .eq. infty) then ncm=nc-1 else ncm=nc endif c c Initialize quantities which will be accumulated c during the SIS iterations. c totvv=0.0d0 totv2=0.0d0 pmax=0.0d0 if (dosrv) then muth=0.0d0 varth=0.0d0 do 3 k=1,ntim surv(k)=0.0d0 sigma(k)=0.0d0 3 continue if(domcse)then do 31 i=1,ntim semu(i)=0.0d0 v2s2(i)=0.0d0 v2s1(i)=0.0d0 31 continue endif else do 333 k=1,ntim do 444 j=1,6 beta(j,k)=0.0d0 444 continue 333 continue mutot=nn+totm0 if (time(1).eq.0.0) then tbeg=2 else tbeg=1 endif c endif c c c Begin SIS iterations. c do 900 jrep=1,nrep c c do 10 i=1,nn amass(i)=0 anxt(i)=0 aval(i)=0.0 10 continue do 20 k=1,ncm cmass(k)=0 cm2(k)=0 cbeg(k)=0 cend(k)=0 20 continue c c Generate first data value. c x=rpari(aa0,bb0,lo(1),hi(1)) c c The predictive probability for the first observation c is NOT random. So we just set it to 1. vv=1.0d0 aval(1)=x aaii=aa0+1.0d0 bbii=bb0+x amass(1)=1 ck=cindex(x,inlo(1),inhi(1),cc,nc) c The next statement is needed only for ci's. cpos(1)=ck cm2(ck)=1 cmass(ck)=1 cbeg(ck)=1 cend(ck)=1 c c Generate other data values. c apos=2 do 200 ii=2,nn iimass=0 do 50 k=inlo(ii),inhi(ii) iimass=iimass+cmass(k) 50 continue mcont=totm0*ppari(aaii,bbii,lo(ii),hi(ii)) vv=vv*(mcont+iimass)/(totm0+ii-1) if (ranf() .le. mcont/(mcont+iimass)) then x=rpari(aaii,bbii,lo(ii),hi(ii)) aval(apos)=x aaii=aaii+1.0 bbii=bbii+x amass(apos)=1 ck=cindex(x,inlo(ii),inhi(ii),cc,nc) c The next statement is needed only for ci's. cpos(apos)=ck cm2(ck)=cm2(ck)+1 if (cmass(ck).eq.0) then cbeg(ck)=apos else anxt(cend(ck))=apos endif cend(ck)=apos apos=apos+1 else u=ranf()*iimass cum=0.0 do 60 k=inlo(ii),inhi(ii) cum=cum+cmass(k) if (u .le. cum) then ck=k goto 61 endif 60 continue 61 continue u=ranf()*cmass(ck) a=cbeg(ck) cum=amass(a) 69 if(u .gt. cum) then a=anxt(a) cum=cum+amass(a) goto 69 endif amass(a)=amass(a)+1 c Accumulate sum of squared atom masses in each interval. cm2(ck)=cm2(ck)+2*amass(a)-1 endif cmass(ck)=cmass(ck)+1 200 continue c c Generate theta. c theta=rgam(aaii,bbii) c c Accumulate survival function info and other output values. c totvv=totvv+vv totv2=totv2+vv*vv pmax=dmax1(pmax,vv) if (dosrv) then muth=muth+vv*theta varth=varth+vv*theta*theta tail=0 c c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c prob(1)=exp(-theta*time(1)) tmpsrv(1)=0.0d0 do 650 i=2,ntim tmpsrv(i)=0.0d0 prob(i)=exp(-theta*time(i)) prob(i-1)=prob(i-1)-prob(i) 650 continue do 652 ck=1,ncm if(cmass(ck).ge.1)then cm=cmass(ck) bot=intim(ck) top=intim(ck+1) if(bot.eq.(top-1))then tmpsrv(bot)=tmpsrv(bot)+cm c Note: sigma does not have to be incremented in this case. else totp=pexpi(theta,time(bot),time(top)) cump=0.0 do 690 i=top-1,bot+1,-1 intp=prob(i)/totp cump=cump+intp sigma(i)=sigma(i)+vv*cm2(ck)*cump*(1-cump) tmpsrv(i)=tmpsrv(i)+cm*intp 690 continue tmpsrv(bot)=tmpsrv(bot)+cm*prob(bot)/totp endif endif 652 continue tmpsrv(ntim)=tmpsrv(ntim)+totm0*prob(ntim) surv(ntim)=surv(ntim)+vv*tmpsrv(ntim) sigma(ntim)=sigma(ntim)+ $ vv*(tmpsrv(ntim)+tmpsrv(ntim)**2) do 730 i=ntim-1,1,-1 tmpsrv(i)=tmpsrv(i)+tmpsrv(i+1)+totm0*prob(i) surv(i)=surv(i)+vv*tmpsrv(i) sigma(i)=sigma(i)+vv*(tmpsrv(i)+tmpsrv(i)**2) 730 continue if(domcse) then do 735 i=1,ntim tmp=vv*tmpsrv(i) v2s1(i)=v2s1(i)+vv*tmp v2s2(i)=v2s2(i)+tmp*tmp 735 continue endif c else c do 6500 i=1,ntim mu(i)=0.0 6500 continue do 7000 ap=1,apos-1 cp=cpos(ap) bot=intim(cp) top=intim(cp+1)-1 if(bot.eq.top)then mu(bot)=mu(bot)+amass(ap) else tp=cindex(aval(ap),bot,top,time,ntim) mu(tp)=mu(tp)+amass(ap) endif 7000 continue do 7300 i=ntim-1,1,-1 mu(i)=mu(i)+mu(i+1) 7300 continue do 731 i=1,ntim mu(i)=mu(i)+totm0*exp(-theta*time(i)) 731 continue c c Now compute the beta cdf, etc. needed to compute the c probability intervals. c do 821 i=tbeg,ntim call vbeta(guess(1,i),mu(i),mutot-mu(i),tol,btmp) do 822 j=1,6 beta(j,i)=beta(j,i)+vv*btmp(j) 822 continue 821 continue c c End the case where doci = .TRUE. endif c c End of SIS iteration. Now start another. c 900 continue c c Now divide through by totvv and normalize various quantities. c pmax=pmax/totvv ess=totvv**2/totv2 if (dosrv) then muth=muth/totvv varth=varth/totvv - muth*muth den=nn+totm0 den2=den*(den+1) do 950 k=1,ntim surv(k)=surv(k)/(totvv*den) sigma(k)=sigma(k)/(totvv*den2) - surv(k)**2 sigma(k)=dsqrt(dmax1(0.0d0,sigma(k))) 950 continue if(domcse)then tmp=den*den tmp2=totvv*totvv do 955 i=1,ntim v2s1(i)=v2s1(i)/den v2s2(i)=v2s2(i)/tmp semu(i)=v2s2(i) + (surv(i)**2)*totv2 $ - 2.0*surv(i)*v2s1(i) semu(i)=dsqrt(dmax1(0.0d0,semu(i)/tmp2)) 955 continue endif else c do 9100 j=1,6 do 911 k=tbeg,ntim beta(j,k)=beta(j,k)/totvv 911 continue 9100 continue c c End case for doci = .TRUE.. c endif c return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function rpari(a,b,t1,t2) implicit none c This function generates a random shifted Pareto with c parameters a,b conditional on it lying in the interval (t1,t2). double precision a,b,t1,t2 double precision infty,tie,tmp real ranf parameter(infty= -99.0,tie= -33.0) c if (t2 .eq. infty) then c xxx Is this next line responsible for loss of accuracy? rpari = t1+(b+t1)*(-1.0+exp(-log(ranf())/a)) else if (t2 .gt. t1) then tmp = 1 - (1 - exp(a*log((b+t1)/(b+t2))))*ranf() rpari = t1+(b+t1)*(-1.0+exp(-log(tmp)/a)) else if (t2 .eq. tie) then rpari = t1 else c write(*,*) "rpari: bad values for t1, t2." rpari=infty endif return end double precision function ppari(a,b,t1,t2) implicit none double precision a,b,t1,t2 double precision infty,tie parameter(infty= -99.0,tie= -33.0) c if (t2 .eq. infty) then ppari = exp(a*log(b/(b+t1))) else if (t2 .gt. t1) then ppari = exp(a*log(b/(b+t1))) - exp(a*log(b/(b+t2))) else if (t2 .eq. tie) then ppari = 0.0 else c write(*,*) "ppari: bad values for t1, t2." ppari=infty endif return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Incorporates modifications to ensure uniform ergodicity. c Also, updates observations in nested groups. c c This subroutine implements the Gibbs sampler based on the Polya c urn scheme for the simplest c situation: a Gamma(a,b) mixture of Dirichlet processes whose c alpha measures are exponential(theta) distributions. The total c mass alpha(R) is assumed constant, that is, the mass does c not depend on the value of theta. c c subroutine gibbs2(lo,hi,inlo,inhi,nn,aa0,bb0,totm0, $ nrep,g0,g1,gextra,cc,intim,nc,time,ntim,new, $ surv,sigma,mu,var) c c Declarations for subroutine arguments. c implicit none integer nn,nc,ntim integer nrep,inlo(nn),inhi(nn),intim(nc),g0,g1,gextra double precision lo(nn),hi(nn),aa0,bb0, $ totm0,cc(nc),time(ntim+1) double precision surv(ntim),sigma(ntim),mu,var logical new c c *** INPUT VALUES *** c lo, hi: Vectors containing the interval information. c nn: the number of data observations. c Observation i lies in the interval (lo(i), hi(i)]. c aa0, bb0: the parameters of the Gamma prior for theta. c totm0: the fixed value of alpha(R). (totm0 is short for total mass.) c nrep: the number of simulation repetitions desired. c new: If new is .true., we assume this is a new gibbs run and c the program makes up a starting state. If new is .false., then c the program continues from the last state it reached. c c *** OUTPUT VALUES *** c time: the distinct endpoints of the intervals in the data c plus the additional time points in adtim. These are c ordered from smallest to largest. c ntim: the length of the array time. c surv: the posterior mean of the survival function (at each c of the time points in the array time). c sigma: the posterior standard deviation of the survival function. c mu, var: the posterior mean and variance of theta. c c c save integer nmax,dmx,mxtim double precision infty,tie,big parameter(nmax=210,dmx=2*nmax,infty= -99.0,tie= -33.0, $ mxtim=3*nmax,big=1.0d100) integer ncm integer amass(nmax+1),anxt(nmax+1),aprv(nmax+1),apos(nmax), $ cmass(dmx),cbeg(dmx),cend(dmx),cpos(nmax+1), $ i,jrep,jj,k,ii,iimass,ap,cp, $ zbeg,zend,g,isect(nmax+1,2),cpold,kk, $ bot,top,gpsiz(nmax),ngp,iipos,igp logical first,anew(nmax+1) double precision aval(nmax+1),mcont,theta,u,cum, $ prob(mxtim),tmpsrv(mxtim),am,am2,intp,cump,totp, $ den,den2,aamix,bbmix,ccmx c c Declaring types for functions used. c integer cindex double precision rexpi,pexpi,rgam real ranf c if (new) then if (cc(nc) .eq. infty) then ncm=nc-1 else ncm=nc endif ccmx=cc(ncm) den=(nn+totm0)*nrep den2=den*(nn+totm0+1) c c Create "fake" complete data to be the starting point c for the Gibbs sampler. c All the data values will be distinct (or at least c treated as if they are distinct) and lie c in the appropriate intervals. c Just take value to be the midpoint of the interval c or add one to the lower endpoint if upper endpoint is c infinite. c c Also divide observations into nested groups. c do 20 k=1,nc cmass(k)=0 20 continue do 33 i=1,nn if (hi(i) .eq. infty) hi(i)=big 33 continue ngp=1 gpsiz(1)=1 do 34 i=2,nn if ((lo(i) .le. lo(i-1)) .and. (hi(i) .ge. hi(i-1))) then gpsiz(ngp)=gpsiz(ngp)+1 else ngp=ngp+1 gpsiz(ngp)=1 endif 34 continue do 35 i=1,nn if (hi(i) .eq. big) hi(i)=infty 35 continue do 111 i=1,nn if (hi(i) .eq. infty) then aval(i)=lo(i)+1.0 else aval(i)=(lo(i)+hi(i))/2.0 endif amass(i)=1 anew(i)=.true. cp=cindex(aval(i),inlo(i),inhi(i),cc,nc) apos(i)=i cpos(i)=cp if(cmass(cp) .eq. 0) then cbeg(cp)=i cend(cp)=i else aprv(i)=cend(cp) anxt(cend(cp))=i cend(cp)=i endif cmass(cp)=cmass(cp)+1 111 continue amass(nn+1)=0 anew(nn+1)=.true. zbeg=nn+1 zend=nn+1 endif c c The following quantities are initialized to zero every c time this subroutine is called, even when new=.false.. c mu=0.0d0 var=0.0d0 do 3 k=1,ntim surv(k)=0.0d0 sigma(k)=0.0d0 3 continue c c c Begin Gibbs iterations (sweeps). c The iterations 1 to g0-1 are discarded. This is the warm-up. c After that the output statistics are recorded (updated) every c g1 iterations. So g0 is the warm-up and g1 is the spacing. c Every gextra iterations of the "basic" gibbs sampler, we do the c "extra step". A single iteration of the "complete" Gibbs c will then involve gextra iterations of the basic Gibbs sampler c followed by an "extra step". After every g1 iterations of the c complete sampler, we record the output statistics. This means c that output statistics are recorded every g1*gextra iterations c of the basic sampler. c first=.true. g=g0 do 900 jrep=1,nrep do 890 jj=1,g do 880 kk=1,gextra c c c Generate theta. c aamix=aa0 bbmix=bb0 do 57 ap=1,nn+1 if (amass(ap).ge.1) then if (cpos(ap).eq.ncm) then bbmix=bbmix+ccmx else aamix=aamix+1.0d0 bbmix=bbmix+aval(ap) endif endif 57 continue theta=rgam(aamix,bbmix) c c Generate the "complete" data values. c iipos=0 do 200 igp=1,ngp c c Delete the observations in this group. c do 198 ii=iipos+1,iipos+gpsiz(igp) ap=apos(ii) cp=cpos(ap) amass(ap)=amass(ap)-1 cmass(cp)=cmass(cp)-1 if(amass(ap) .eq. 0) then if (cmass(cp) .ge. 1) then if (cbeg(cp) .eq. ap) then cbeg(cp)=anxt(ap) else if (cend(cp) .eq. ap) then cend(cp)=aprv(ap) else anxt(aprv(ap))=anxt(ap) aprv(anxt(ap))=aprv(ap) endif endif anxt(zend)=ap zend=ap endif 198 continue c c Generate new values for this group. c do 199 ii=iipos+1,iipos+gpsiz(igp) iimass=0 do 50 k=inlo(ii),inhi(ii) iimass=iimass+cmass(k) 50 continue mcont=totm0*pexpi(theta,lo(ii),hi(ii)) if (ranf() .le. mcont/(mcont+iimass)) then ap=zbeg zbeg=anxt(zbeg) aval(ap)=rexpi(theta,lo(ii),hi(ii)) amass(ap)=1 cp=cindex(aval(ap),inlo(ii),inhi(ii),cc,nc) cpos(ap)=cp if (cmass(cp).eq.0) then cbeg(cp)=ap else anxt(cend(cp))=ap aprv(ap)=cend(cp) endif cend(cp)=ap else u=ranf()*iimass cum=0.0 do 60 k=inlo(ii),inhi(ii) cum=cum+cmass(k) if (u .le. cum) then cp=k goto 61 endif 60 continue 61 continue u=ranf()*cmass(cp) ap=cbeg(cp) cum=amass(ap) 69 if (u .gt. cum) then ap=anxt(ap) cum=cum+amass(ap) goto 69 endif amass(ap)=amass(ap)+1 endif apos(ii)=ap cmass(cp)=cmass(cp)+1 199 continue iipos=iipos+gpsiz(igp) 200 continue c 880 continue c c Now add the "extra step". Re-randomize the positions c of the atoms, respecting the current group structure. c Before doing this, we generate a new value of theta in a way c designed to ensure uniform ergodicity. c do 2900 ii=1,nn ap=apos(ii) if (anew(ap)) then isect(ap,1)=inlo(ii) isect(ap,2)=inhi(ii) anew(ap)=.false. else isect(ap,1)=max(isect(ap,1),inlo(ii)) isect(ap,2)=min(isect(ap,2),inhi(ii)) endif 2900 continue aamix=aa0 bbmix=bb0 do 2910 ap=1,nn+1 if (amass(ap).ge.1) then anew(ap)=.true. if (cc(isect(ap,2)+1).eq.infty) then bbmix=bbmix+cc(isect(ap,1)) else aamix=aamix+1.0d0 bbmix=bbmix+aval(ap) endif endif 2910 continue theta=rgam(aamix,bbmix) do 2920 ap=1,nn+1 if (amass(ap).ge.1) then aval(ap)=rexpi(theta,cc(isect(ap,1)),cc(isect(ap,2)+1)) if(isect(ap,2).gt.isect(ap,1))then cpold=cpos(ap) cp=cindex(aval(ap),isect(ap,1),isect(ap,2),cc,nc) if(cp.ne.cpold)then cpos(ap)=cp cmass(cpold)=cmass(cpold)-amass(ap) if(cmass(cpold) .ge. 1) then if(cbeg(cpold) .eq. ap) then cbeg(cpold)=anxt(ap) else if(cend(cpold) .eq. ap) then cend(cpold)=aprv(ap) else anxt(aprv(ap))=anxt(ap) aprv(anxt(ap))=aprv(ap) endif endif if (cmass(cp).eq.0) then cbeg(cp)=ap else anxt(cend(cp))=ap aprv(ap)=cend(cp) endif cend(cp)=ap cmass(cp)=cmass(cp)+amass(ap) endif endif endif 2920 continue c 890 continue c if (first .eqv. .true.) then g=g1 first=.false. endif c c Accumulate survival function info and other output values. c This uses "Fancy" Rao-Blackwellization which involves c conditioning on the group structure. c mu=mu+theta var=var+theta*theta c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c prob(1)=exp(-theta*time(1)) tmpsrv(1)=0.0 do 650 i=2,ntim tmpsrv(i)=0.0 prob(i)=exp(-theta*time(i)) prob(i-1)=prob(i-1)-prob(i) 650 continue do 700 ap=1,nn+1 if (amass(ap).ge.1) then bot=intim(isect(ap,1)) top=intim(isect(ap,2)+1) if(bot.eq.(top-1))then tmpsrv(bot)=tmpsrv(bot)+amass(ap) c Note: sigma does not have to be incremented in this case. else totp=pexpi(theta,time(bot),time(top)) cump=0.0 am=amass(ap) am2=am*am do 690 i=top-1,bot+1,-1 intp=prob(i)/totp cump=cump+intp sigma(i)=sigma(i)+am2*cump*(1-cump) tmpsrv(i)=tmpsrv(i)+am*intp 690 continue tmpsrv(bot)=tmpsrv(bot)+am*prob(bot)/totp endif endif 700 continue tmpsrv(ntim)=tmpsrv(ntim)+totm0*prob(ntim) surv(ntim)=surv(ntim)+tmpsrv(ntim) sigma(ntim)=sigma(ntim)+tmpsrv(ntim)+tmpsrv(ntim)**2 do 730 i=ntim-1,1,-1 tmpsrv(i)=tmpsrv(i)+tmpsrv(i+1)+totm0*prob(i) surv(i)=surv(i)+tmpsrv(i) sigma(i)=sigma(i)+tmpsrv(i)+tmpsrv(i)**2 730 continue c c 900 continue c c Now convert totals to means (normalize various quantities). c do 910 k=1,ntim surv(k)=surv(k)/den sigma(k)=sigma(k)/den2 c c Warning: The above line does not actually compute the posterior sigma. c We must convert this into the true sigma (in the Splus program). To c compute sigma in this program, use the following two lines of code: c sigma(k)=sigma(k)/den2-surv(k)**2 c sigma(k)=dsqrt(dmax1(0.0d0,sigma(k))) c 910 continue mu=mu/nrep var=var/nrep c Warning: This must be converted to var by subtracting mu*mu as c done below. c var=var/nrep - mu*mu c c call intpr("batch done",-1,0,1) return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine grnest(lo,hi,nn,ngp) c c Arrange the observations so the data intervals c are nested in groups. c c Declarations for subroutine arguments. c implicit none integer nn,ngp double precision lo(nn),hi(nn) c integer nmax,tn,i,j double precision infty,blank,big,glo,ghi parameter(nmax=210,infty=-99.0,blank=-888.0,big=1.0d100) double precision tlo(nmax),thi(nmax) c do 3 i=1,nn if (hi(i) .eq. infty) hi(i)=big 3 continue c tn=0 ngp=0 5 if (tn .lt. nn) then do 10 i=1,nn if (lo(i).ne.blank) then tn=tn+1 ngp=ngp+1 glo=lo(i) ghi=hi(i) tlo(tn)=glo thi(tn)=ghi lo(i)=blank j=i goto 11 endif 10 continue 11 continue do 20 i=j+1,nn if (lo(i).ne.blank) then if ((lo(i).le.glo).and.(hi(i).ge.ghi)) then tn=tn+1 glo=lo(i) ghi=hi(i) tlo(tn)=glo thi(tn)=ghi lo(i)=blank endif endif 20 continue goto 5 endif c do 30 i=1,nn lo(i)=tlo(i) hi(i)=thi(i) 30 continue do 35 i=1,nn if (hi(i) .eq. big) hi(i)=infty 35 continue return end cccccccccccccccccccccccccccccc subroutine vbeta(y,a,b,tol,v) c c Compute the cdf, pdf and the derivative of the pdf c for the beta(a,b) distribution at the value y. c implicit none double precision y(2),a,b,tol,v(6) c real w,w1,betaln double precision bln,yy,tolc integer ierr c tolc=1.0d0-tol bln=betaln(real(a),real(b)) c yy=y(1) if ((yy .gt. tol).and.(yy .lt. tolc)) then call bratio(real(a),real(b),real(yy),real(1.0d0-yy), $ w,w1,ierr) v(1)=w v(2)=exp(-bln+(a-1.0)*log(yy)+(b-1.0)*log(1.0d0-yy)) v(3)=((a-1.0)/yy - (b-1.0)/(1.0-yy))*v(2) endif c yy=y(2) if ((yy .gt. tol).and.(yy .lt. tolc)) then call bratio(real(b),real(a),real(yy),real(1.0d0-yy), $ w,w1,ierr) v(4)=w v(5)=exp(-bln+(b-1.0)*log(yy)+(a-1.0)*log(1.0d0-yy)) v(6)=((b-1.0)/yy - (a-1.0)/(1.0-yy))*v(5) endif return end ccccccccccccccc subroutine impgss(yy,nn,p,tol,ans) c y[1,] is the guess for the confidence limit. c y[2,] is the estimated lower tail probability at y[1,]. c y[3,] is the estimated density at y[1,]. c y[4,] is the estimate derivative at y[1,]. c Assume log H(y) = a*log(y) + b + c*y. Solve for a, b, c. c This is a naive Newton-Raphson scheme for solving c log(y)+ d*y = e. implicit none integer nn double precision yy(4,nn),p,tol,ans(nn) integer i,nsh,niter,mxit double precision aa,bb,cc,eps,y,d,e,fun,der,step,oy, $ ofun,diff,odiff,s1,newy,ub,lb,tolc parameter(mxit=20) parameter(eps=1.0d-8) tolc=1.0d0-tol do 30 i=1,nn y=yy(1,i) if ((y .gt. tol).and.(y .lt. tolc)) then aa = y**2 *((yy(3,i)/yy(2,i))**2-yy(4,i)/yy(2,i)) cc = yy(3,i)/yy(2,i) - aa/y bb = log(yy(2,i)) - aa*log(y) - cc*y c Solve a*log(y) + b + c*y = log(p) c or equivalently log(y) + (c/a)*y = (log(p)-b)/a. d = cc/aa e = (log(p)-bb)/aa ub=min(2.0*y,(1.0+y)/2.0) lb=y/2.0 c If possible, make the initial guess is to the left of the c global maximum when d < 0 and aa>0. c If aa<0 the function is behaving strangely, and it is c probably best just to go to the nearest root. if ((aa.gt.0).and.(d.lt.0.0).and.(y.ge.-1.0/d)) then y=-0.5/d if (y.gt.ub) then y=ub else if (y.lt.lb) then y=lb endif endif fun=log(y)+d*y diff=abs(fun-e) niter=0 10 if((diff.ge.eps).and.(niter.lt.mxit))then niter=niter+1 oy=y ofun=fun odiff=diff der=(1/oy)+d step=(e-ofun)/der newy=oy+step if(((y.lt.lb+eps).and.(newy.lt.lb)).or. $ ((y.gt.ub-eps).and.(newy.gt.ub))) goto 19 if (newy.lt.lb) then step=lb-oy else if (newy.gt.ub) then step=ub-oy endif if(niter.eq.1)s1=step nsh=0 12 if(nsh.lt.10)then y=oy+step fun=log(y)+d*y diff=abs(fun-e) c Make sure absolute difference actually decreases. if (diff.ge.odiff)then step=step/2 nsh=nsh+1 goto 12 endif endif goto 10 endif 19 if (niter.ge.mxit) then ans(i)=yy(1,i)+0.5*s1 else ans(i)=y endif else ans(i)=yy(1,i) endif 30 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Parametric Bayes model. (Limit as totm goes to infinity.) c c subroutine parmod(lo,hi,nn,aa0,bb0,nrep,g0,g1,gextra, $ time,ntim,new,surv,sigma,mu,var) c c Declarations for subroutine arguments. c implicit none integer nn,ntim integer nrep,g0,g1,gextra double precision lo(nn),hi(nn),aa0,bb0,time(ntim+1) double precision surv(ntim),sigma(ntim),mu,var logical new c c *** INPUT VALUES *** c lo, hi: Vectors containing the interval information. c nn: the number of data observations. c Observation i lies in the interval (lo(i), hi(i)]. c aa0, bb0: the parameters of the Gamma prior for theta. c nrep: the number of simulation repetitions desired. c new: If new is .true., we assume this is a new gibbs run and c the program makes up a starting state. If new is .false., then c the program continues from the last state it reached. c c *** OUTPUT VALUES *** c time: the distinct endpoints of the intervals in the data c plus the additional time points in adtim. These are c ordered from smallest to largest. c ntim: the length of the array time. c surv: the posterior mean of the survival function (at each c of the time points in the array time). c sigma: the posterior standard deviation of the survival function. c mu, var: the posterior mean and variance of theta. c c c save integer nmax,dmx,resumd,mxtim double precision infty,tie parameter(nmax=210,dmx=2*nmax,infty= -99.0,tie= -33.0, $ resumd=40,mxtim=3*nmax) integer i,jrep,jj,k,ii,nsumd,g,kk logical first double precision aval(nmax+1),theta, $ sumdis,aa,bb c c Declaring types for functions used. c double precision rexpi,rgam c if(new)then aa=aa0+nn nsumd=0 c c Create "fake" complete data to be the starting point c for the Gibbs sampler. c All the data values will be distinct (or at least c treated as if they are distinct) and lie c in the appropriate intervals. c Just take value to be the midpoint of the interval c or add one to the lower endpoint if upper endpoint is c infinite. c sumdis=0.0d0 do 111 i=1,nn if (hi(i) .eq. infty) then aval(i)=lo(i)+1.0 else aval(i)=(lo(i)+hi(i))/2.0 endif sumdis=sumdis+aval(i) 111 continue endif c c The following quantities are initialized to zero every c time this subroutine is called, even when new=.false.. c mu=0.0d0 var=0.0d0 do 3 k=1,ntim surv(k)=0.0d0 sigma(k)=0.0d0 3 continue c c c Begin Gibbs iterations (sweeps). c The iterations 1 to g0-1 are discarded. This is the warm-up. c After that the output statistics are recorded (updated) every c g1 iterations. So g0 is the warm-up and g1 is the spacing. c Every gextra iterations of the "basic" gibbs sampler, we do the c "extra step". A single iteration of the "complete" Gibbs c will then involve gextra iterations of the basic Gibbs sampler c followed by an "extra step". After every g1 iterations of the c complete sampler, we record the output statistics. This means c that output statistics are recorded every g1*gextra iterations c of the basic sampler. c first=.true. g=g0 do 900 jrep=1,nrep do 890 jj=1,g do 880 kk=1,gextra c c c Generate theta. c theta=rgam(aa,bb0+sumdis) c c Generate the "complete" data values. c do 200 ii=1,nn sumdis=sumdis-aval(ii) aval(ii)=rexpi(theta,lo(ii),hi(ii)) sumdis=sumdis+aval(ii) 200 continue c c Every resumd iterations recompute the value of sumdis. c This guards against accumulation of round-off error. c nsumd=nsumd+1 if (nsumd .eq. resumd) then nsumd=0 sumdis=0.0d0 do 747 ii=1,nn sumdis=sumdis+aval(ii) 747 continue endif c 880 continue c 890 continue c if (first .eqv. .true.) then g=g1 first=.false. endif c c mu=mu+theta var=var+theta*theta c c This chunk of code updates the posterior variance c (as well as the posterior mean) of the survival function. c bb=bb0+sumdis do 650 i=1,ntim surv(i)=surv(i)+exp(aa*log(bb/(bb+time(i)))) sigma(i)=sigma(i)+exp(aa*log(bb/(bb+2.0*time(i)))) 650 continue c 900 continue c c Now convert totals to means (normalize various quantities). c do 910 k=1,ntim surv(k)=surv(k)/nrep sigma(k)=sigma(k)/nrep c c Warning: The above line does not actually compute the posterior sigma. c We must convert this into the true sigma (in the Splus program). To c compute sigma in this program, use the following two lines of code: c sigma(k)=sigma(k)/den2-surv(k)**2 c sigma(k)=dsqrt(dmax1(0.0d0,sigma(k))) c 910 continue mu=mu/nrep var=var/nrep c Warning: This must be converted to var by subtracting mu*mu as c done below. c var=var/nrep - mu*mu c return end subroutine vranf(n,v) implicit none integer n,i double precision v(n) real ranf c do 10 i=1,n v(i)=ranf() 10 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccc subroutines written by others cccccccccccccccccc ccccccccccccccccc (sometimes slightly modified by me) ccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C ALGORITHM 654, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 13, NO. 3, P. 318. SUBROUTINE GRATIO (A, X, ANS, QANS, IND) C ---------------------------------------------------------------------- C EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS C P(A,X) AND Q(A,X) C C ---------- C C IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X C ARE NOT BOTH 0. C C ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE C P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER. C IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS C POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF C IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE C 6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY C IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT. C C ERROR RETURN ... C ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE, C WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT. C P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN C X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE. C ---------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C -------------------- REAL J, L, ACC0(3), BIG(3), E00(3), X00(3), WK(20) REAL D0(13), D1(12), D2(10), D3(8), D4(6), D5(4), D6(2) C -------------------- DATA ACC0(1)/5.E-15/, ACC0(2)/5.E-7/, ACC0(3)/5.E-4/ DATA BIG(1)/20.0/, BIG(2)/14.0/, BIG(3)/10.0/ DATA E00(1)/.25E-3/, E00(2)/.25E-1/, E00(3)/.14/ DATA X00(1)/31.0/, X00(2)/17.0/, X00(3)/9.7/ C -------------------- C ALOG10 = LN(10) C RT2PIN = 1/SQRT(2*PI) C RTPI = SQRT(PI) C -------------------- DATA ALOG10/2.30258509299405/ DATA RT2PIN/.398942280401433/ DATA RTPI /1.77245385090552/ DATA THIRD /.333333333333333/ C -------------------- DATA D0(1) / .833333333333333E-01/, D0(2) /-.148148148148148E-01/, * D0(3) / .115740740740741E-02/, D0(4) / .352733686067019E-03/, * D0(5) /-.178755144032922E-03/, D0(6) / .391926317852244E-04/, * D0(7) /-.218544851067999E-05/, D0(8) /-.185406221071516E-05/, * D0(9) / .829671134095309E-06/, D0(10)/-.176659527368261E-06/, * D0(11)/ .670785354340150E-08/, D0(12)/ .102618097842403E-07/, * D0(13)/-.438203601845335E-08/ C -------------------- DATA D10 /-.185185185185185E-02/, D1(1) /-.347222222222222E-02/, * D1(2) / .264550264550265E-02/, D1(3) /-.990226337448560E-03/, * D1(4) / .205761316872428E-03/, D1(5) /-.401877572016461E-06/, * D1(6) /-.180985503344900E-04/, D1(7) / .764916091608111E-05/, * D1(8) /-.161209008945634E-05/, D1(9) / .464712780280743E-08/, * D1(10)/ .137863344691572E-06/, D1(11)/-.575254560351770E-07/, * D1(12)/ .119516285997781E-07/ C -------------------- DATA D20 / .413359788359788E-02/, D2(1) /-.268132716049383E-02/, * D2(2) / .771604938271605E-03/, D2(3) / .200938786008230E-05/, * D2(4) /-.107366532263652E-03/, D2(5) / .529234488291201E-04/, * D2(6) /-.127606351886187E-04/, D2(7) / .342357873409614E-07/, * D2(8) / .137219573090629E-05/, D2(9) /-.629899213838006E-06/, * D2(10)/ .142806142060642E-06/ C -------------------- DATA D30 / .649434156378601E-03/, D3(1) / .229472093621399E-03/, * D3(2) /-.469189494395256E-03/, D3(3) / .267720632062839E-03/, * D3(4) /-.756180167188398E-04/, D3(5) /-.239650511386730E-06/, * D3(6) / .110826541153473E-04/, D3(7) /-.567495282699160E-05/, * D3(8) / .142309007324359E-05/ C -------------------- DATA D40 /-.861888290916712E-03/, D4(1) / .784039221720067E-03/, * D4(2) /-.299072480303190E-03/, D4(3) /-.146384525788434E-05/, * D4(4) / .664149821546512E-04/, D4(5) /-.396836504717943E-04/, * D4(6) / .113757269706784E-04/ C -------------------- DATA D50 /-.336798553366358E-03/, D5(1) /-.697281375836586E-04/, * D5(2) / .277275324495939E-03/, D5(3) /-.199325705161888E-03/, * D5(4) / .679778047793721E-04/ C -------------------- DATA D60 / .531307936463992E-03/, D6(1) /-.592166437353694E-03/, * D6(2) / .270878209671804E-03/ C -------------------- DATA D70 / .344367606892378E-03/ C -------------------- C ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST C FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 . C E = SPMPAR(1) C C -------------------- IF (A .LT. 0.0 .OR. X .LT. 0.0) GO TO 400 IF (A .EQ. 0.0 .AND. X .EQ. 0.0) GO TO 400 IF (A*X .EQ. 0.0) GO TO 331 C IOP = IND + 1 IF (IOP .NE. 1 .AND. IOP .NE. 2) IOP = 3 ACC = AMAX1(ACC0(IOP),E) E0 = E00(IOP) X0 = X00(IOP) C C SELECT THE APPROPRIATE ALGORITHM C IF (A .GE. 1.0) GO TO 10 IF (A .EQ. 0.5) GO TO 320 IF (X .LT. 1.1) GO TO 110 T1 = A*ALOG(X) - X U = A*EXP(T1) IF (U .EQ. 0.0) GO TO 310 R = U*(1.0 + GAM1(A)) GO TO 170 C 10 IF (A .GE. BIG(IOP)) GO TO 20 IF (A .GT. X .OR. X .GE. X0) GO TO 11 TWOA = A + A M = INT(TWOA) IF (TWOA .NE. FLOAT(M)) GO TO 11 I = M/2 IF (A .EQ. FLOAT(I)) GO TO 140 GO TO 150 11 T1 = A*ALOG(X) - X R = EXP(T1)/GAMMA(A) GO TO 30 C 20 L = X/A IF (L .EQ. 0.0) GO TO 300 S = 0.5 + (0.5 - L) Z = RLOG(L) IF (Z .GE. 700.0/A) GO TO 330 Y = A*Z RTA = SQRT(A) IF (ABS(S) .LE. E0/RTA) GO TO 250 IF (ABS(S) .LE. 0.4) GO TO 200 C T = (1.0/A)**2 T1 = (((0.75*T - 1.0)*T + 3.5)*T - 105.0)/(A*1260.0) T1 = T1 - Y R = RT2PIN*RTA*EXP(T1) C 30 IF (R .EQ. 0.0) GO TO 331 IF (X .LE. AMAX1(A,ALOG10)) GO TO 50 IF (X .LT. X0) GO TO 170 GO TO 80 C C TAYLOR SERIES FOR P/R C 50 APN = A + 1.0 T = X/APN WK(1) = T DO 51 N = 2,20 APN = APN + 1.0 T = T*(X/APN) IF (T .LE. 1.E-3) GO TO 60 51 WK(N) = T N = 20 C 60 SUM = T TOL = 0.5*ACC 61 APN = APN + 1.0 T = T*(X/APN) SUM = SUM + T IF (T .GT. TOL) GO TO 61 C MAX = N - 1 DO 70 M = 1,MAX N = N - 1 70 SUM = SUM + WK(N) ANS = (R/A)*(1.0 + SUM) QANS = 0.5 + (0.5 - ANS) RETURN C C ASYMPTOTIC EXPANSION C 80 AMN = A - 1.0 T = AMN/X WK(1) = T DO 81 N = 2,20 AMN = AMN - 1.0 T = T*(AMN/X) IF (ABS(T) .LE. 1.E-3) GO TO 90 81 WK(N) = T N = 20 C 90 SUM = T 91 IF (ABS(T) .LE. ACC) GO TO 100 AMN = AMN - 1.0 T = T*(AMN/X) SUM = SUM + T GO TO 91 C 100 MAX = N - 1 DO 101 M = 1,MAX N = N - 1 101 SUM = SUM + WK(N) QANS = (R/X)*(1.0 + SUM) ANS = 0.5 + (0.5 - QANS) RETURN C C TAYLOR SERIES FOR P(A,X)/X**A C 110 AN = 3.0 C = X SUM = X/(A + 3.0) TOL = 3.0*ACC/(A + 1.0) 111 AN = AN + 1.0 C = -C*(X/AN) T = C/(A + AN) SUM = SUM + T IF (ABS(T) .GT. TOL) GO TO 111 J = A*X*((SUM/6.0 - 0.5/(A + 2.0))*X + 1.0/(A + 1.0)) C Z = A*ALOG(X) H = GAM1(A) G = 1.0 + H IF (X .LT. 0.25) GO TO 120 IF (A .LT. X/2.59) GO TO 135 GO TO 130 120 IF (Z .GT. -.13394) GO TO 135 C 130 W = EXP(Z) ANS = W*G*(0.5 + (0.5 - J)) QANS = 0.5 + (0.5 - ANS) RETURN C 135 L = REXP(Z) W = 0.5 + (0.5 + L) QANS = (W*J - L)*G - H IF (QANS .LT. 0.0) GO TO 310 ANS = 0.5 + (0.5 - QANS) RETURN C C FINITE SUMS FOR Q WHEN A .GE. 1 C AND 2*A IS AN INTEGER C 140 SUM = EXP(-X) T = SUM N = 1 C = 0.0 GO TO 160 C 150 RTX = SQRT(X) SUM = ERFC1(0,RTX) T = EXP(-X)/(RTPI*RTX) N = 0 C = -0.5 C 160 IF (N .EQ. I) GO TO 161 N = N + 1 C = C + 1.0 T = (X*T)/C SUM = SUM + T GO TO 160 161 QANS = SUM ANS = 0.5 + (0.5 - QANS) RETURN C C CONTINUED FRACTION EXPANSION C 170 TOL = AMAX1(5.0*E,ACC) A2NM1 = 1.0 A2N = 1.0 B2NM1 = X B2N = X + (1.0 - A) C = 1.0 171 A2NM1 = X*A2N + C*A2NM1 B2NM1 = X*B2N + C*B2NM1 AM0 = A2NM1/B2NM1 C = C + 1.0 CMA = C - A A2N = A2NM1 + CMA*A2N B2N = B2NM1 + CMA*B2N AN0 = A2N/B2N IF (ABS(AN0 - AM0) .GE. TOL*AN0) GO TO 171 C QANS = R*AN0 ANS = 0.5 + (0.5 - QANS) RETURN C C GENERAL TEMME EXPANSION C 200 IF (ABS(S) .LE. 2.0*E .AND. A*E*E .GT. 3.28E-3) GO TO 400 C = EXP(-Y) W = 0.5*ERFC1(1,SQRT(Y)) U = 1.0/A Z = SQRT(Z + Z) IF (L .LT. 1.0) Z = -Z IF (IOP - 2) 210,220,230 C 210 IF (ABS(S) .LE. 1.E-3) GO TO 260 C0 = ((((((((((((D0(13) * Z + D0(12)) * Z + D0(11)) * Z * + D0(10)) * Z + D0(9)) * Z + D0(8)) * Z + D0(7)) * Z * + D0(6)) * Z + D0(5)) * Z + D0(4)) * Z + D0(3)) * Z * + D0(2)) * Z + D0(1)) * Z - THIRD C1 = (((((((((((D1(12) * Z + D1(11)) * Z + D1(10)) * Z * + D1(9)) * Z + D1(8)) * Z + D1(7)) * Z + D1(6)) * Z * + D1(5)) * Z + D1(4)) * Z + D1(3)) * Z + D1(2)) * Z * + D1(1)) * Z + D10 C2 = (((((((((D2(10) * Z + D2(9)) * Z + D2(8)) * Z * + D2(7)) * Z + D2(6)) * Z + D2(5)) * Z + D2(4)) * Z * + D2(3)) * Z + D2(2)) * Z + D2(1)) * Z + D20 C3 = (((((((D3(8) * Z + D3(7)) * Z + D3(6)) * Z * + D3(5)) * Z + D3(4)) * Z + D3(3)) * Z + D3(2)) * Z * + D3(1)) * Z + D30 C4 = (((((D4(6) * Z + D4(5)) * Z + D4(4)) * Z + D4(3)) * Z * + D4(2)) * Z + D4(1)) * Z + D40 C5 = (((D5(4) * Z + D5(3)) * Z + D5(2)) * Z + D5(1)) * Z * + D50 C6 = (D6(2) * Z + D6(1)) * Z + D60 T = ((((((D70*U + C6)*U + C5)*U + C4)*U + C3)*U + C2)*U * + C1)*U + C0 GO TO 240 C 220 C0 = (((((D0(6) * Z + D0(5)) * Z + D0(4)) * Z + D0(3)) * Z * + D0(2)) * Z + D0(1)) * Z - THIRD C1 = (((D1(4) * Z + D1(3)) * Z + D1(2)) * Z + D1(1)) * Z * + D10 C2 = D2(1) * Z + D20 T = (C2*U + C1)*U + C0 GO TO 240 C 230 T = ((D0(3) * Z + D0(2)) * Z + D0(1)) * Z - THIRD C 240 IF (L .LT. 1.0) GO TO 241 QANS = C*(W + RT2PIN*T/RTA) ANS = 0.5 + (0.5 - QANS) RETURN 241 ANS = C*(W - RT2PIN*T/RTA) QANS = 0.5 + (0.5 - ANS) RETURN C C TEMME EXPANSION FOR L = 1 C 250 IF (A*E*E .GT. 3.28E-3) GO TO 400 C = 0.5 + (0.5 - Y) W = (0.5 - SQRT(Y)*(0.5 + (0.5 - Y/3.0))/RTPI)/C U = 1.0/A Z = SQRT(Z + Z) IF (L .LT. 1.0) Z = -Z IF (IOP - 2) 260,270,280 C 260 C0 = ((((((D0(7) * Z + D0(6)) * Z + D0(5)) * Z + D0(4)) * Z * + D0(3)) * Z + D0(2)) * Z + D0(1)) * Z - THIRD C1 = (((((D1(6) * Z + D1(5)) * Z + D1(4)) * Z + D1(3)) * Z * + D1(2)) * Z + D1(1)) * Z + D10 C2 = ((((D2(5) * Z + D2(4)) * Z + D2(3)) * Z + D2(2)) * Z * + D2(1)) * Z + D20 C3 = (((D3(4) * Z + D3(3)) * Z + D3(2)) * Z + D3(1)) * Z * + D30 C4 = (D4(2) * Z + D4(1)) * Z + D40 C5 = (D5(2) * Z + D5(1)) * Z + D50 C6 = D6(1) * Z + D60 T = ((((((D70*U + C6)*U + C5)*U + C4)*U + C3)*U + C2)*U * + C1)*U + C0 GO TO 240 C 270 C0 = (D0(2) * Z + D0(1)) * Z - THIRD C1 = D1(1) * Z + D10 T = (D20*U + C1)*U + C0 GO TO 240 C 280 T = D0(1) * Z - THIRD GO TO 240 C C SPECIAL CASES C 300 ANS = 0.0 QANS = 1.0 RETURN C 310 ANS = 1.0 QANS = 0.0 RETURN C 320 IF (X .GE. 0.25) GO TO 321 ANS = ERFff(SQRT(X)) QANS = 0.5 + (0.5 - ANS) RETURN 321 QANS = ERFC1(0,SQRT(X)) ANS = 0.5 + (0.5 - QANS) RETURN C 330 IF (ABS(S) .LE. 2.0*E) GO TO 400 331 IF (X .LE. A) GO TO 300 GO TO 310 C C ERROR RETURN C 400 ANS = 2.0 RETURN END SUBROUTINE GAMINV (A, X, X0, P, Q, IERR) C ---------------------------------------------------------------------- C INVERSE INCOMPLETE GAMMA RATIO FUNCTION C C GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1. C THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER C ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X C TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE C PARTICULAR COMPUTER ARITHMETIC BEING USED. C C ------------ C C X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0, C AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT C NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN C A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE C IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X. C C X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER C DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET C X0 .LE. 0. C C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. C WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING C VALUES ... C C IERR = 0 THE SOLUTION WAS OBTAINED. ITERATION WAS C NOT USED. C IERR.GT.0 THE SOLUTION WAS OBTAINED. IERR ITERATIONS C WERE PERFORMED. C IERR = -2 (INPUT ERROR) A .LE. 0 C IERR = -3 NO SOLUTION WAS OBTAINED. THE RATIO Q/A C IS TOO LARGE. C IERR = -4 (INPUT ERROR) P + Q .NE. 1 C IERR = -6 20 ITERATIONS WERE PERFORMED. THE MOST C RECENT VALUE OBTAINED FOR X IS GIVEN. C THIS CANNOT OCCUR IF X0 .LE. 0. C IERR = -7 ITERATION FAILED. NO VALUE IS GIVEN FOR X. C THIS MAY OCCUR WHEN X IS APPROXIMATELY 0. C IERR = -8 A VALUE FOR X HAS BEEN OBTAINED, BUT THE C ROUTINE IS NOT CERTAIN OF ITS ACCURACY. C ITERATION CANNOT BE PERFORMED IN THIS C CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY C WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS C POSITIVE THEN THIS CAN OCCUR WHEN A IS C EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY C LARGE (SAY A .GE. 1.E20). C ---------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C ------------------- REAL LN10, EPS0(2), AMIN(2), BMIN(2), DMIN(2), EMIN(2) C ------------------- C LN10 = LN(10) C C = EULER CONSTANT C ------------------- DATA LN10 /2.302585/ DATA C /.577215664901533/ C ------------------- DATA A0 /3.31125922108741/, A1 /11.6616720288968/, * A2 /4.28342155967104/, A3 /.213623493715853/ DATA B1 /6.61053765625462/, B2 /6.40691597760039/, * B3 /1.27364489782223/, B4 /.036117081018842/ C ------------------- DATA EPS0(1) /1.E-10/, EPS0(2) /1.E-08/ DATA AMIN(1) / 500.0/, AMIN(2) / 100.0/ DATA BMIN(1) /1.E-28/, BMIN(2) /1.E-13/ DATA DMIN(1) /1.E-06/, DMIN(2) /1.E-04/ DATA EMIN(1) /2.E-03/, EMIN(2) /6.E-03/ C ------------------- DATA TOL /1.E-5/ C ------------------- C ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS. C E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0. C XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE C LARGEST POSITIVE NUMBER. C E = SPMPAR(1) XMIN = SPMPAR(2) XMAX = SPMPAR(3) C ------------------- X = 0.0 IF (A .LE. 0.0) GO TO 500 T = DBLE(P) + DBLE(Q) - 1.D0 IF (ABS(T) .GT. E) GO TO 520 C IERR = 0 IF (P .EQ. 0.0) RETURN IF (Q .EQ. 0.0) GO TO 400 IF (A .EQ. 1.0) GO TO 410 C E2 = 2.0*E AMAX = 0.4E-10/(E*E) IOP = 1 IF (E .GT. 1.E-10) IOP = 2 EPS = EPS0(IOP) XN = X0 IF (X0 .GT. 0.0) GO TO 100 C C SELECTION OF THE INITIAL APPROXIMATION XN OF X C WHEN A .LT. 1 C IF (A .GT. 1.0) GO TO 30 G = GAMMA(A + 1.0) QG = Q*G IF (QG .EQ. 0.0) GO TO 560 B = QG/A IF (QG .GT. 0.6*A) GO TO 20 IF (A .GE. 0.30 .OR. B .LT. 0.35) GO TO 10 T = EXP(-(B + C)) U = T*EXP(T) XN = T*EXP(U) GO TO 100 C 10 IF (B .GE. 0.45) GO TO 20 IF (B .EQ. 0.0) GO TO 560 Y = -ALOG(B) S = 0.5 + (0.5 - A) Z = ALOG(Y) T = Y - S*Z IF (B .LT. 0.15) GO TO 11 XN = Y - S*ALOG(T) - ALOG(1.0 + S/(T + 1.0)) GO TO 200 11 IF (B .LE. 0.01) GO TO 12 U = ((T + 2.0*(3.0 - A))*T + (2.0 - A)*(3.0 - A))/ * ((T + (5.0 - A))*T + 2.0) XN = Y - S*ALOG(T) - ALOG(U) GO TO 200 12 C1 = -S*Z C2 = -S*(1.0 + C1) C3 = S*((0.5*C1 + (2.0 - A))*C1 + (2.5 - 1.5*A)) C4 = -S*(((C1/3.0 + (2.5 - 1.5*A))*C1 + ((A - 6.0)*A + 7.0))*C1 * + ((11.0*A - 46)*A + 47.0)/6.0) C5 = -S*((((-C1/4.0 + (11.0*A - 17.0)/6.0)*C1 * + ((-3.0*A + 13.0)*A - 13.0))*C1 * + 0.5*(((2.0*A - 25.0)*A + 72.0)*A - 61.0))*C1 * + (((25.0*A - 195.0)*A + 477.0)*A -379.0)/12.0) XN = ((((C5/Y + C4)/Y + C3)/Y + C2)/Y + C1) + Y IF (A .GT. 1.0) GO TO 200 IF (B .GT. BMIN(IOP)) GO TO 200 X = XN RETURN C 20 IF (B*Q .GT. 1.E-8) GO TO 21 XN = EXP(-(Q/A + C)) GO TO 23 21 IF (P .LE. 0.9) GO TO 22 XN = EXP((ALNREL(-Q) + GAMLN1(A))/A) GO TO 23 22 XN = EXP(ALOG(P*G)/A) 23 IF (XN .EQ. 0.0) GO TO 510 T = 0.5 + (0.5 - XN/(A + 1.0)) XN = XN/T GO TO 100 C C SELECTION OF THE INITIAL APPROXIMATION XN OF X C WHEN A .GT. 1 C 30 IF (Q .LE. 0.5) GO TO 31 W = ALOG(P) GO TO 32 31 W = ALOG(Q) 32 T = SQRT(-2.0*W) S = T - (((A3*T + A2)*T + A1)*T + A0)/((((B4*T + B3)*T * + B2)*T + B1)*T + 1.0) IF (Q .GT. 0.5) S = -S C RTA = SQRT(A) S2 = S*S XN = A + S*RTA + (S2 - 1.0)/3.0 + S*(S2 - 7.0)/(36.0*RTA) 1 - ((3.0*S2 + 7.0)*S2 - 16.0)/(810.0*A) 2 + S*((9.0*S2 + 256.0)*S2 - 433.0)/(38880.0*A*RTA) XN = AMAX1(XN, 0.0) IF (A .LT. AMIN(IOP)) GO TO 40 X = XN D = 0.5 + (0.5 - X/A) IF (ABS(D) .LE. DMIN(IOP)) RETURN C 40 IF (P .LE. 0.5) GO TO 50 IF (XN .LT. 3.0*A) GO TO 200 Y = -(W + GAMLN(A)) D = AMAX1(2.0, A*(A - 1.0)) IF (Y .LT. LN10*D) GO TO 41 S = 1.0 - A Z = ALOG(Y) GO TO 12 41 T = A - 1.0 XN = Y + T*ALOG(XN) - ALNREL(-T/(XN + 1.0)) XN = Y + T*ALOG(XN) - ALNREL(-T/(XN + 1.0)) GO TO 200 C 50 AP1 = A + 1.0 IF (XN .GT. 0.70*AP1) GO TO 101 W = W + GAMLN(AP1) IF (XN .GT. 0.15*AP1) GO TO 60 AP2 = A + 2.0 AP3 = A + 3.0 X = EXP((W + X)/A) X = EXP((W + X - ALOG(1.0 + (X/AP1)*(1.0 + X/AP2)))/A) X = EXP((W + X - ALOG(1.0 + (X/AP1)*(1.0 + X/AP2)))/A) X = EXP((W + X - ALOG(1.0 + (X/AP1)*(1.0 + (X/AP2)*(1.0 * + X/AP3))))/A) XN = X IF (XN .GT. 1.E-2*AP1) GO TO 60 IF (XN .LE. EMIN(IOP)*AP1) RETURN GO TO 101 C 60 APN = AP1 T = XN/APN SUM = 1.0 + T 61 APN = APN + 1.0 T = T*(XN/APN) SUM = SUM + T IF (T .GT. 1.E-4) GO TO 61 T = W - ALOG(SUM) XN = EXP((XN + T)/A) XN = XN*(1.0 - (A*ALOG(XN) - XN -T)/(A - XN)) GO TO 101 C C SCHRODER ITERATION USING P C 100 IF (P .GT. 0.5) GO TO 200 101 IF (P .LE. 1.E10*XMIN) GO TO 550 AM1 = (A - 0.5) - 0.5 102 IF (A .LE. AMAX) GO TO 110 D = 0.5 + (0.5 - XN/A) IF (ABS(D) .LE. E2) GO TO 550 C 110 IF (IERR .GE. 20) GO TO 530 IERR = IERR + 1 CALL GRATIO(A,XN,PN,QN,0) IF (PN .EQ. 0.0 .OR. QN .EQ. 0.0) GO TO 550 R = RCOMP(A,XN) IF (R .EQ. 0.0) GO TO 550 T = (PN - P)/R W = 0.5*(AM1 - XN) IF (ABS(T) .LE. 0.1 .AND. ABS(W*T) .LE. 0.1) GO TO 120 X = XN*(1.0 - T) IF (X .LE. 0.0) GO TO 540 D = ABS(T) GO TO 121 C 120 H = T*(1.0 + W*T) X = XN*(1.0 - H) IF (X .LE. 0.0) GO TO 540 IF (ABS(W) .GE. 1.0 .AND. ABS(W)*T*T .LE. EPS) RETURN D = ABS(H) 121 XN = X IF (D .GT. TOL) GO TO 102 IF (D .LE. EPS) RETURN IF (ABS(P - PN) .LE. TOL*P) RETURN GO TO 102 C C SCHRODER ITERATION USING Q C 200 IF (Q .LE. 1.E10*XMIN) GO TO 550 AM1 = (A - 0.5) - 0.5 201 IF (A .LE. AMAX) GO TO 210 D = 0.5 + (0.5 - XN/A) IF (ABS(D) .LE. E2) GO TO 550 C 210 IF (IERR .GE. 20) GO TO 530 IERR = IERR + 1 CALL GRATIO(A,XN,PN,QN,0) IF (PN .EQ. 0.0 .OR. QN .EQ. 0.0) GO TO 550 R = RCOMP(A,XN) IF (R .EQ. 0.0) GO TO 550 T = (Q - QN)/R W = 0.5*(AM1 - XN) IF (ABS(T) .LE. 0.1 .AND. ABS(W*T) .LE. 0.1) GO TO 220 X = XN*(1.0 - T) IF (X .LE. 0.0) GO TO 540 D = ABS(T) GO TO 221 C 220 H = T*(1.0 + W*T) X = XN*(1.0 - H) IF (X .LE. 0.0) GO TO 540 IF (ABS(W) .GE. 1.0 .AND. ABS(W)*T*T .LE. EPS) RETURN D = ABS(H) 221 XN = X IF (D .GT. TOL) GO TO 201 IF (D .LE. EPS) RETURN IF (ABS(Q - QN) .LE. TOL*Q) RETURN GO TO 201 C C SPECIAL CASES C 400 X = XMAX RETURN C 410 IF (Q .LT. 0.9) GO TO 411 X = -ALNREL(-P) RETURN 411 X = -ALOG(Q) RETURN C C ERROR RETURN C 500 IERR = -2 RETURN C 510 IERR = -3 RETURN C 520 IERR = -4 RETURN C 530 IERR = -6 RETURN C 540 IERR = -7 RETURN C 550 X = XN IERR = -8 RETURN C 560 X = XMAX IERR = -8 RETURN END FUNCTION ERFff(X) C ****************************************************************** C EVALUATION OF THE REAL ERROR FUNCTION C ****************************************************************** DIMENSION A(4),B(4),P(8),Q(8),R(5),S(5) DATA A(1)/-1.65581836870402E-4/, A(2)/3.25324098357738E-2/, 1 A(3)/1.02201136918406E-1/, A(4)/1.12837916709552E00/ DATA B(1)/4.64988945913179E-3/, B(2)/7.01333417158511E-2/, 1 B(3)/4.23906732683201E-1/, B(4)/1.00000000000000E00/ DATA P(1)/-1.36864857382717E-7/, P(2)/5.64195517478974E-1/, 1 P(3)/7.21175825088309E00/, P(4)/4.31622272220567E01/, 2 P(5)/1.52989285046940E02/, P(6)/3.39320816734344E02/, 3 P(7)/4.51918953711873E02/, P(8)/3.00459261020162E02/ DATA Q(1)/1.00000000000000E00/, Q(2)/1.27827273196294E01/, 1 Q(3)/7.70001529352295E01/, Q(4)/2.77585444743988E02/, 2 Q(5)/6.38980264465631E02/, Q(6)/9.31354094850610E02/, 3 Q(7)/7.90950925327898E02/, Q(8)/3.00459260956983E02/ DATA R(1)/2.10144126479064E00/, R(2)/2.62370141675169E01/, 1 R(3)/2.13688200555087E01/, R(4)/4.65807828718470E00/, 2 R(5)/2.82094791773523E-1/ DATA S(1)/9.41537750555460E01/, S(2)/1.87114811799590E02/, 1 S(3)/9.90191814623914E01/, S(4)/1.80124575948747E01/, 2 S(5)/1.00000000000000E00/ DATA C/5.64189583547756E-1/ C ------------------- AX=ABS(X) X2=AX*AX IF (AX.GE.0.5) GO TO 10 TOP=((A(1)*X2+A(2))*X2+A(3))*X2+A(4) BOT=((B(1)*X2+B(2))*X2+B(3))*X2+B(4) ERFff=X*TOP/BOT RETURN C 10 IF (AX.GT.4.0) GO TO 20 TOP=((((((P(1)*AX+P(2))*AX+P(3))*AX+P(4))*AX+P(5))*AX * +P(6))*AX+P(7))*AX+P(8) BOT=((((((Q(1)*AX+Q(2))*AX+Q(3))*AX+Q(4))*AX+Q(5))*AX * +Q(6))*AX+Q(7))*AX+Q(8) ERFff=1.0-EXP(-X2)*TOP/BOT IF (X.LT.0.0) ERFff=-ERFff RETURN C 20 ERFff=1.0 IF (AX.GE.5.54) GO TO 21 T=1.0/X2 TOP=(((R(1)*T+R(2))*T+R(3))*T+R(4))*T+R(5) BOT=(((S(1)*T+S(2))*T+S(3))*T+S(4))*T+S(5) ERFff=C-TOP/(X2*BOT) ERFff=1.0-EXP(-X2)*ERFff/AX 21 IF (X.LT.0.0) ERFff=-ERFff RETURN END REAL FUNCTION ERFC1(IND,X) C ---------------------------------------------------------------------- C EVALUATION OF THE REAL COMPLEMENTARY ERROR FUNCTION C C ERFC1(IND,X) = ERFC(X) IF IND = 0 C ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE C ---------------------------------------------------------------------- DIMENSION A(4),B(4),P(8),Q(8),R(5),S(5) DATA A(1)/-1.65581836870402E-4/, A(2)/3.25324098357738E-2/, 1 A(3)/1.02201136918406E-1/, A(4)/1.12837916709552E00/ DATA B(1)/4.64988945913179E-3/, B(2)/7.01333417158511E-2/, 1 B(3)/4.23906732683201E-1/, B(4)/1.00000000000000E00/ DATA P(1)/-1.36864857382717E-7/, P(2)/5.64195517478974E-1/, 1 P(3)/7.21175825088309E00/, P(4)/4.31622272220567E01/, 2 P(5)/1.52989285046940E02/, P(6)/3.39320816734344E02/, 3 P(7)/4.51918953711873E02/, P(8)/3.00459261020162E02/ DATA Q(1)/1.00000000000000E00/, Q(2)/1.27827273196294E01/, 1 Q(3)/7.70001529352295E01/, Q(4)/2.77585444743988E02/, 2 Q(5)/6.38980264465631E02/, Q(6)/9.31354094850610E02/, 3 Q(7)/7.90950925327898E02/, Q(8)/3.00459260956983E02/ DATA R(1)/2.10144126479064E00/, R(2)/2.62370141675169E01/, 1 R(3)/2.13688200555087E01/, R(4)/4.65807828718470E00/, 2 R(5)/2.82094791773523E-1/ DATA S(1)/9.41537750555460E01/, S(2)/1.87114811799590E02/, 1 S(3)/9.90191814623914E01/, S(4)/1.80124575948747E01/, 2 S(5)/1.00000000000000E00/ DATA C/5.64189583547756E-1/ C ------------------- AX=ABS(X) X2=AX*AX IF (AX.GE.0.47) GO TO 10 TOP=((A(1)*X2+A(2))*X2+A(3))*X2+A(4) BOT=((B(1)*X2+B(2))*X2+B(3))*X2+B(4) ERFC1=1.0-X*TOP/BOT IF (IND.NE.0) ERFC1=EXP(X2)*ERFC1 RETURN C 10 IF (AX.GT.4.0) GO TO 20 TOP=((((((P(1)*AX+P(2))*AX+P(3))*AX+P(4))*AX+P(5))*AX * +P(6))*AX+P(7))*AX+P(8) BOT=((((((Q(1)*AX+Q(2))*AX+Q(3))*AX+Q(4))*AX+Q(5))*AX * +Q(6))*AX+Q(7))*AX+Q(8) ERFC1=TOP/BOT IF (IND.EQ.0) GO TO 11 IF (X.LT.0.0) ERFC1=2.0*EXP(X2)-ERFC1 RETURN 11 ERFC1=EXP(-X2)*ERFC1 IF (X.LT.0.0) ERFC1=2.0-ERFC1 RETURN C 20 IF (X.LE.-5.33) GO TO 30 T=1.0/X2 TOP=(((R(1)*T+R(2))*T+R(3))*T+R(4))*T+R(5) BOT=(((S(1)*T+S(2))*T+S(3))*T+S(4))*T+S(5) ERFC1=(C-TOP/(X2*BOT))/AX IF (IND.EQ.0) GO TO 11 IF (X.LT.0.0) ERFC1=2.0*EXP(X2)-ERFC1 RETURN C 30 ERFC1=2.0 IF (IND.NE.0) ERFC1=EXP(X2)*ERFC1 RETURN END REAL FUNCTION REXP(X) C ------------------------------------------------------------------ C COMPUTATION OF EXP(X) - 1 C ------------------------------------------------------------------ DATA P1/ .914041914819518E-09/, P2/ .238082361044469E-01/, * Q1/-.499999999085958E+00/, Q2/ .107141568980644E+00/, * Q3/-.119041179760821E-01/, Q4/ .595130811860248E-03/ C ------------------ IF (ABS(X) .GT. 0.15) GO TO 10 REXP = X*(((P2*X + P1)*X + 1.0)/((((Q4*X + Q3)*X + Q2)*X * + Q1)*X + 1.0)) RETURN C 10 W = EXP(X) IF (X .GT. 0.0) GO TO 20 REXP = (W - 0.5) - 0.5 RETURN 20 REXP = W*(0.5 + (0.5 - 1.0/W)) RETURN END REAL FUNCTION ALNREL(A) C ------------------------------------------------------------------ C EVALUATION OF THE FUNCTION LN(1 + A) C ------------------------------------------------------------------ DATA P1/-.129418923021993E+01/, P2/.405303492862024E+00/, * P3/-.178874546012214E-01/ DATA Q1/-.162752256355323E+01/, Q2/.747811014037616E+00/, * Q3/-.845104217945565E-01/ C --------------------- IF (ABS(A) .GT. 0.375) GO TO 10 T = A/(A + 2.0) T2 = T*T W = (((P3*T2 + P2)*T2 + P1)*T2 + 1.0)/ * (((Q3*T2 + Q2)*T2 + Q1)*T2 + 1.0) ALNREL = 2.0*T*W RETURN C 10 X = 1.D0 + DBLE(A) ALNREL = ALOG(X) RETURN END REAL FUNCTION RLOG(X) C ------------------- C COMPUTATION OF X - 1 - LN(X) C ------------------- DATA A/.566749439387324E-01/ DATA B/.456512608815524E-01/ C ------------------- DATA P0/ .333333333333333E+00/, P1/-.224696413112536E+00/, * P2/ .620886815375787E-02/ DATA Q1/-.127408923933623E+01/, Q2/ .354508718369557E+00/ C ------------------- IF (X .LT. 0.61 .OR. X .GT. 1.57) GO TO 100 IF (X .LT. 0.82) GO TO 10 IF (X .GT. 1.18) GO TO 20 C C ARGUMENT REDUCTION C U = (X - 0.5) - 0.5 W1 = 0.0 GO TO 30 C 10 U = DBLE(X) - 0.7D0 U = U/0.7 W1 = A - U*0.3 GO TO 30 C 20 U = 0.75D0*DBLE(X) - 1.D0 W1 = B + U/3.0 C C SERIES EXPANSION C 30 R = U/(U + 2.0) T = R*R W = ((P2*T + P1)*T + P0)/((Q2*T + Q1)*T + 1.0) RLOG = 2.0*T*(1.0/(1.0 - R) - R*W) + W1 RETURN C C 100 R = (X - 0.5) - 0.5 RLOG = R - ALOG(X) RETURN END REAL FUNCTION GAMMA(A) C----------------------------------------------------------------------- C C EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS C C ----------- C C GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT C BE COMPUTED. C C----------------------------------------------------------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C----------------------------------------------------------------------- REAL P(7), Q(7) DOUBLE PRECISION D, G, Z, LNX, GLOG C-------------------------- C D = 0.5*(LN(2*PI) - 1) C-------------------------- DATA PI /3.1415926535898/ DATA D /.41893853320467274178D0/ C-------------------------- DATA P(1)/ .539637273585445E-03/, P(2)/ .261939260042690E-02/, 1 P(3)/ .204493667594920E-01/, P(4)/ .730981088720487E-01/, 2 P(5)/ .279648642639792E+00/, P(6)/ .553413866010467E+00/, 3 P(7)/ 1.0/ DATA Q(1)/-.832979206704073E-03/, Q(2)/ .470059485860584E-02/, 1 Q(3)/ .225211131035340E-01/, Q(4)/-.170458969313360E+00/, 2 Q(5)/-.567902761974940E-01/, Q(6)/ .113062953091122E+01/, 3 Q(7)/ 1.0/ C-------------------------- DATA R1/.820756370353826E-03/, R2/-.595156336428591E-03/, 1 R3/.793650663183693E-03/, R4/-.277777777770481E-02/, 2 R5/.833333333333333E-01/ C-------------------------- GAMMA = 0.0 X = A IF (ABS(A) .GE. 15.0) GO TO 60 C----------------------------------------------------------------------- C EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15 C----------------------------------------------------------------------- T = 1.0 M = INT(A) - 1 C C LET T BE THE PRODUCT OF A-J WHEN A .GE. 2 C IF (M) 20,12,10 10 DO 11 J = 1,M X = X - 1.0 11 T = X*T 12 X = X - 1.0 GO TO 40 C C LET T BE THE PRODUCT OF A+J WHEN A .LT. 1 C 20 T = A IF (A .GT. 0.0) GO TO 30 M = - M - 1 IF (M .EQ. 0) GO TO 22 DO 21 J = 1,M X = X + 1.0 21 T = X*T 22 X = (X + 0.5) + 0.5 T = X*T IF (T .EQ. 0.0) RETURN C 30 CONTINUE C C THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS C CODE MAY BE OMITTED IF DESIRED. C IF (ABS(T) .GE. 1.E-30) GO TO 40 IF (ABS(T)*SPMPAR(3) .LE. 1.0001) RETURN GAMMA = 1.0/T RETURN C C COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1 C 40 TOP = P(1) BOT = Q(1) DO 41 I = 2,7 TOP = P(I) + X*TOP 41 BOT = Q(I) + X*BOT GAMMA = TOP/BOT C C TERMINATION C IF (A .LT. 1.0) GO TO 50 GAMMA = GAMMA*T RETURN 50 GAMMA = GAMMA/T RETURN C----------------------------------------------------------------------- C EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15 C----------------------------------------------------------------------- 60 IF (ABS(A) .GE. 1.E3) RETURN IF (A .GT. 0.0) GO TO 70 X = -A N = X T = X - N IF (T .GT. 0.9) T = 1.0 - T S = SIN(PI*T)/PI IF (MOD(N,2) .EQ. 0) S = -S IF (S .EQ. 0.0) RETURN C C COMPUTE THE MODIFIED ASYMPTOTIC SUM C 70 T = 1.0/(X*X) G = ((((R1*T + R2)*T + R3)*T + R4)*T + R5)/X C C ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X) C BUT LESS ACCURACY WILL NORMALLY BE OBTAINED. C LNX = GLOG(X) C C FINAL ASSEMBLY C Z = X G = (D + G) + (Z - 0.5D0)*(LNX - 1.D0) W = G T = G - DBLE(W) IF (W .GT. 0.99999*EXPARG(0)) RETURN GAMMA = EXP(W)*(1.0 + T) IF (A .LT. 0.0) GAMMA = (1.0/(GAMMA*S))/X RETURN END DOUBLE PRECISION FUNCTION GLOG(X) C ------------------- C EVALUATION OF LN(X) FOR X .GE. 15 C ------------------- REAL X DOUBLE PRECISION Z, W(163) C ------------------- DATA C1/.286228750476730/, C2/.399999628131494/, 1 C3/.666666666752663/ C ------------------- C W(J) = LN(J + 14) FOR EACH J C ------------------- DATA W(1) /.270805020110221007D+01/, 1 W(2) /.277258872223978124D+01/, W(3) /.283321334405621608D+01/, 2 W(4) /.289037175789616469D+01/, W(5) /.294443897916644046D+01/, 3 W(6) /.299573227355399099D+01/, W(7) /.304452243772342300D+01/, 4 W(8) /.309104245335831585D+01/, W(9) /.313549421592914969D+01/, 5 W(10)/.317805383034794562D+01/, W(11)/.321887582486820075D+01/, 6 W(12)/.325809653802148205D+01/, W(13)/.329583686600432907D+01/, 7 W(14)/.333220451017520392D+01/, W(15)/.336729582998647403D+01/, 8 W(16)/.340119738166215538D+01/, W(17)/.343398720448514625D+01/, 9 W(18)/.346573590279972655D+01/, W(19)/.349650756146648024D+01/, 1 W(20)/.352636052461616139D+01/, W(21)/.355534806148941368D+01/, 2 W(22)/.358351893845611000D+01/, W(23)/.361091791264422444D+01/, 3 W(24)/.363758615972638577D+01/, W(25)/.366356164612964643D+01/, 4 W(26)/.368887945411393630D+01/, W(27)/.371357206670430780D+01/, 5 W(28)/.373766961828336831D+01/, W(29)/.376120011569356242D+01/, 6 W(30)/.378418963391826116D+01/ DATA W(31)/.380666248977031976D+01/, 1 W(32)/.382864139648909500D+01/, W(33)/.385014760171005859D+01/, 2 W(34)/.387120101090789093D+01/, W(35)/.389182029811062661D+01/, 3 W(36)/.391202300542814606D+01/, W(37)/.393182563272432577D+01/, 4 W(38)/.395124371858142735D+01/, W(39)/.397029191355212183D+01/, 5 W(40)/.398898404656427438D+01/, W(41)/.400733318523247092D+01/, 6 W(42)/.402535169073514923D+01/, W(43)/.404305126783455015D+01/, 7 W(44)/.406044301054641934D+01/, W(45)/.407753744390571945D+01/, 8 W(46)/.409434456222210068D+01/, W(47)/.411087386417331125D+01/, 9 W(48)/.412713438504509156D+01/, W(49)/.414313472639153269D+01/, 1 W(50)/.415888308335967186D+01/, W(51)/.417438726989563711D+01/, 2 W(52)/.418965474202642554D+01/, W(53)/.420469261939096606D+01/, 3 W(54)/.421950770517610670D+01/, W(55)/.423410650459725938D+01/, 4 W(56)/.424849524204935899D+01/, W(57)/.426267987704131542D+01/, 5 W(58)/.427666611901605531D+01/, W(59)/.429045944114839113D+01/, 6 W(60)/.430406509320416975D+01/ DATA W(61)/.431748811353631044D+01/, 1 W(62)/.433073334028633108D+01/, W(63)/.434380542185368385D+01/, 2 W(64)/.435670882668959174D+01/, W(65)/.436944785246702149D+01/, 3 W(66)/.438202663467388161D+01/, W(67)/.439444915467243877D+01/, 4 W(68)/.440671924726425311D+01/, W(69)/.441884060779659792D+01/, 5 W(70)/.443081679884331362D+01/, W(71)/.444265125649031645D+01/, 6 W(72)/.445434729625350773D+01/, W(73)/.446590811865458372D+01/, 7 W(74)/.447733681447820647D+01/, W(75)/.448863636973213984D+01/, 8 W(76)/.449980967033026507D+01/, W(77)/.451085950651685004D+01/, 9 W(78)/.452178857704904031D+01/, W(79)/.453259949315325594D+01/, 1 W(80)/.454329478227000390D+01/, W(81)/.455387689160054083D+01/, 2 W(82)/.456434819146783624D+01/, W(83)/.457471097850338282D+01/, 3 W(84)/.458496747867057192D+01/, W(85)/.459511985013458993D+01/, 4 W(86)/.460517018598809137D+01/, W(87)/.461512051684125945D+01/, 5 W(88)/.462497281328427108D+01/, W(89)/.463472898822963577D+01/, 6 W(90)/.464439089914137266D+01/ DATA W(91) /.465396035015752337D+01/, 1 W(92) /.466343909411206714D+01/, W(93) /.467282883446190617D+01/, 2 W(94) /.468213122712421969D+01/, W(95) /.469134788222914370D+01/, 3 W(96) /.470048036579241623D+01/, W(97) /.470953020131233414D+01/, 4 W(98) /.471849887129509454D+01/, W(99) /.472738781871234057D+01/, 5 W(100)/.473619844839449546D+01/, W(101)/.474493212836325007D+01/, 6 W(102)/.475359019110636465D+01/, W(103)/.476217393479775612D+01/, 7 W(104)/.477068462446566476D+01/, W(105)/.477912349311152939D+01/, 8 W(106)/.478749174278204599D+01/, W(107)/.479579054559674109D+01/, 9 W(108)/.480402104473325656D+01/, W(109)/.481218435537241750D+01/, 1 W(110)/.482028156560503686D+01/, W(111)/.482831373730230112D+01/, 2 W(112)/.483628190695147800D+01/, W(113)/.484418708645859127D+01/, 3 W(114)/.485203026391961717D+01/, W(115)/.485981240436167211D+01/, 4 W(116)/.486753445045558242D+01/, W(117)/.487519732320115154D+01/, 5 W(118)/.488280192258637085D+01/, W(119)/.489034912822175377D+01/, 6 W(120)/.489783979995091137D+01/ DATA W(121)/.490527477843842945D+01/, 1 W(122)/.491265488573605201D+01/, W(123)/.491998092582812492D+01/, 2 W(124)/.492725368515720469D+01/, W(125)/.493447393313069176D+01/, 3 W(126)/.494164242260930430D+01/, W(127)/.494875989037816828D+01/, 4 W(128)/.495582705760126073D+01/, W(129)/.496284463025990728D+01/, 5 W(130)/.496981329957600062D+01/, W(131)/.497673374242057440D+01/, 6 W(132)/.498360662170833644D+01/, W(133)/.499043258677873630D+01/, 7 W(134)/.499721227376411506D+01/, W(135)/.500394630594545914D+01/, 8 W(136)/.501063529409625575D+01/, W(137)/.501727983681492433D+01/, 9 W(138)/.502388052084627639D+01/, W(139)/.503043792139243546D+01/, 1 W(140)/.503695260241362916D+01/, W(141)/.504342511691924662D+01/, 2 W(142)/.504985600724953705D+01/, W(143)/.505624580534830806D+01/, 3 W(144)/.506259503302696680D+01/, W(145)/.506890420222023153D+01/, 4 W(146)/.507517381523382692D+01/, W(147)/.508140436498446300D+01/, 5 W(148)/.508759633523238407D+01/, W(149)/.509375020080676233D+01/, 6 W(150)/.509986642782419842D+01/ DATA W(151)/.510594547390058061D+01/, 1 W(152)/.511198778835654323D+01/, W(153)/.511799381241675511D+01/, 2 W(154)/.512396397940325892D+01/, W(155)/.512989871492307347D+01/, 3 W(156)/.513579843705026176D+01/, W(157)/.514166355650265984D+01/, 4 W(158)/.514749447681345304D+01/, W(159)/.515329159449777895D+01/, 5 W(160)/.515905529921452903D+01/, W(161)/.516478597392351405D+01/, 6 W(162)/.517048399503815178D+01/, W(163)/.517614973257382914D+01/ C IF (X .GE. 178.0) GO TO 10 N = X T = (X - N)/(X + N) T2 = T*T Z = (((C1*T2 + C2)*T2 + C3)*T2 + 2.0)*T GLOG = W(N - 14) + Z RETURN C 10 GLOG = ALOG(X) RETURN END REAL FUNCTION EXPARG (L) C-------------------------------------------------------------------- C COMPUTATION OF THE LARGEST ARGUMENT W FOR WHICH EXP(W) C MAY BE COMPUTED. (ONLY AN APPROXIMATE VALUE IS NEEDED.) C-------------------------------------------------------------------- if (L.eq.0) then EXPARG = 0.99999*ALOG(SPMPAR(3)) else exparg = 0.99999*alog(spmpar(2)) endif RETURN END REAL FUNCTION GAM1(A) C ------------------------------------------------------------------ C COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5 C ------------------------------------------------------------------ REAL P(7), Q(5), R(9) C ------------------- DATA P(1)/ .577215664901533E+00/, P(2)/-.409078193005776E+00/, * P(3)/-.230975380857675E+00/, P(4)/ .597275330452234E-01/, * P(5)/ .766968181649490E-02/, P(6)/-.514889771323592E-02/, * P(7)/ .589597428611429E-03/ C ------------------- DATA Q(1)/ .100000000000000E+01/, Q(2)/ .427569613095214E+00/, * Q(3)/ .158451672430138E+00/, Q(4)/ .261132021441447E-01/, * Q(5)/ .423244297896961E-02/ C ------------------- DATA R(1)/-.422784335098468E+00/, R(2)/-.771330383816272E+00/, * R(3)/-.244757765222226E+00/, R(4)/ .118378989872749E+00/, * R(5)/ .930357293360349E-03/, R(6)/-.118290993445146E-01/, * R(7)/ .223047661158249E-02/, R(8)/ .266505979058923E-03/, * R(9)/-.132674909766242E-03/ C ------------------- DATA S1 / .273076135303957E+00/, S2 / .559398236957378E-01/ C ------------------- T = A D = A - 0.5 IF (D .GT. 0.0) T = D - 0.5 IF (T) 30,10,20 C 10 GAM1 = 0.0 RETURN C 20 TOP = (((((P(7)*T + P(6))*T + P(5))*T + P(4))*T + P(3))*T * + P(2))*T + P(1) BOT = (((Q(5)*T + Q(4))*T + Q(3))*T + Q(2))*T + 1.0 W = TOP/BOT IF (D .GT. 0.0) GO TO 21 GAM1 = A*W RETURN 21 GAM1 = (T/A)*((W - 0.5) - 0.5) RETURN C 30 TOP = (((((((R(9)*T + R(8))*T + R(7))*T + R(6))*T + R(5))*T * + R(4))*T + R(3))*T + R(2))*T + R(1) BOT = (S2*T + S1)*T + 1.0 W = TOP/BOT IF (D .GT. 0.0) GO TO 31 GAM1 = A*((W + 0.5) + 0.5) RETURN 31 GAM1 = T*W/A RETURN END REAL FUNCTION GAMLN(A) C ****************************************************************** C EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A C ****************************************************************** C WRITTEN BY ALFRED H. MORRIS C NAVAL SURFACE WEAPONS CENTER C DAHLGREN, VIRGINIA C --------------------- C D = 0.5*(LN(2*PI) - 1) C --------------------- DATA D/.418938533204673/ C --------------------- DATA C0/.833333333333333E-01/, C1/-.277777777770481E-02/, 1 C2/.793650663183693E-03/, C3/-.595156336428591E-03/, 2 C4/.820756370353826E-03/ C ------------------------------------------------------------------ IF (A .GT. 0.8) GO TO 10 GAMLN = GAMLN1(A) - ALOG(A) RETURN 10 IF (A .GT. 2.25) GO TO 20 X = DBLE(A) - 1.D0 GAMLN = GAMLN1(X) RETURN C 20 IF (A .GE. 15.0) GO TO 30 N = A - 1.25 X = A W = 1.0 DO 21 I = 1,N X = X - 1.0 21 W = X*W GAMLN = GAMLN1(X - 1.0) + ALOG(W) RETURN C 30 X = (1.0/A)**2 W = ((((C4*X + C3)*X + C2)*X + C1)*X + C0)/A GAMLN = (D + W) + (A - 0.5)*(ALOG(A) - 1.0) END REAL FUNCTION GAMLN1(A) C ------------------------------------------------------------------ C EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25 C ------------------------------------------------------------------ DATA P0/ .577215664901533E+00/, P1/ .844203922187225E+00/, * P2/-.168860593646662E+00/, P3/-.780427615533591E+00/, * P4/-.402055799310489E+00/, P5/-.673562214325671E-01/, * P6/-.271935708322958E-02/ DATA Q1/ .288743195473681E+01/, Q2/ .312755088914843E+01/, * Q3/ .156875193295039E+01/, Q4/ .361951990101499E+00/, * Q5/ .325038868253937E-01/, Q6/ .667465618796164E-03/ C ----------------- DATA R0/.422784335098467E+00/, R1/.848044614534529E+00/, * R2/.565221050691933E+00/, R3/.156513060486551E+00/, * R4/.170502484022650E-01/, R5/.497958207639485E-03/ DATA S1/.124313399877507E+01/, S2/.548042109832463E+00/, * S3/.101552187439830E+00/, S4/.713309612391000E-02/, * S5/.116165475989616E-03/ C ----------------- IF (A .GE. 0.6) GO TO 10 W = ((((((P6*A + P5)*A + P4)*A + P3)*A + P2)*A + P1)*A + P0)/ * ((((((Q6*A + Q5)*A + Q4)*A + Q3)*A + Q2)*A + Q1)*A + 1.0) GAMLN1 = -A*W RETURN C 10 X = DBLE(A) - 1.D0 W = (((((R5*X + R4)*X + R3)*X + R2)*X + R1)*X + R0)/ * (((((S5*X + S4)*X + S3)*X + S2)*X + S1)*X + 1.0) GAMLN1 = X*W RETURN END REAL FUNCTION RCOMP(A,X) C ------------------- C EVALUATION OF EXP(-X)*X**A/GAMMA(A) C ------------------- C RT2PIN = 1/SQRT(2*PI) C ------------------- DATA RT2PIN/.398942280401433/ C ------------------- RCOMP = 0.0 IF (A .GE. 20.0) GO TO 20 T = A*ALOG(X) - X IF (A .GE. 1.0) GO TO 10 RCOMP = (A*EXP(T))*(1.0 + GAM1(A)) RETURN 10 RCOMP = EXP(T)/GAMMA(A) RETURN C 20 U = X/A IF (U .EQ. 0.0) RETURN T = (1.0/A)**2 T1 = (((0.75*T - 1.0)*T + 3.5)*T - 105.0)/(A*1260.0) T1 = T1 - A*RLOG(U) RCOMP = RT2PIN*SQRT(A)*EXP(T1) RETURN END real function spmpar(i) integer i c ********** c c Function spmpar c c This function provides single precision machine parameters c when the appropriate set of data statements is activated (by c removing the c from column 1) and all other data statements are c rendered inactive. Most of the parameter values were obtained c from the corresponding Bell Laboratories Port Library function. c c The function statement is c c real function spmpar(i) c c where c c i is an integer input variable set to 1, 2, or 3 which c selects the desired machine parameter. If the machine has c t base b digits and its smallest and largest exponents are c emin and emax, respectively, then these parameters are c c spmpar(1) = b**(1 - t), the machine precision, c c spmpar(2) = b**(emin - 1), the smallest magnitude, c c spmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude. c c Argonne National Laboratory. MINPACK Project. November 1996. c Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More' c c ********** integer mcheps(2) integer minmag(2) integer maxmag(2) real rmach(3) equivalence (rmach(1),mcheps(1)) equivalence (rmach(2),minmag(1)) equivalence (rmach(3),maxmag(1)) c c Machine constants for the IBM 360/370 series, c the Amdahl 470/V6, the ICL 2900, the Itel AS/6, c the Xerox Sigma 5/7/9 and the Sel systems 85/86. c c data rmach(1) / z3c100000 / c data rmach(2) / z00100000 / c data rmach(3) / z7fffffff / c c Machine constants for the Honeywell 600/6000 series. c c data rmach(1) / o716400000000 / c data rmach(2) / o402400000000 / c data rmach(3) / o376777777777 / c c Machine constants for the CDC 6000/7000 series. c c data rmach(1) / 16414000000000000000b / c data rmach(2) / 00014000000000000000b / c data rmach(3) / 37767777777777777777b / c c Machine constants for the PDP-10 (KA or KI processor). c c data rmach(1) / "147400000000 / c data rmach(2) / "000400000000 / c data rmach(3) / "377777777777 / c c Machine constants for the PDP-11 fortran supporting c 32-bit integers (expressed in integer and octal). c c data mcheps(1) / 889192448 / c data minmag(1) / 8388608 / c data maxmag(1) / 2147483647 / c c data rmach(1) / o06500000000 / c data rmach(2) / o00040000000 / c data rmach(3) / o17777777777 / c c Machine constants for the PDP-11 fortran supporting c 16-bit integers (expressed in integer and octal). c c data mcheps(1),mcheps(2) / 13568, 0 / c data minmag(1),minmag(2) / 128, 0 / c data maxmag(1),maxmag(2) / 32767, -1 / c c data mcheps(1),mcheps(2) / o032400, o000000 / c data minmag(1),minmag(2) / o000200, o000000 / c data maxmag(1),maxmag(2) / o077777, o177777 / c c Machine constants for the Burroughs 5700/6700/7700 systems. c c data rmach(1) / o1301000000000000 / c data rmach(2) / o1771000000000000 / c data rmach(3) / o0777777777777777 / c c Machine constants for the Burroughs 1700 system. c c data rmach(1) / z4ea800000 / c data rmach(2) / z400800000 / c data rmach(3) / z5ffffffff / c c Machine constants for the Univac 1100 series. c c data rmach(1) / o147400000000 / c data rmach(2) / o000400000000 / c data rmach(3) / o377777777777 / c c Machine constants for the Data General Eclipse S/200. c c Note - it may be appropriate to include the following card - c static rmach(3) c c data minmag/20k,0/,maxmag/77777k,177777k/ c data mcheps/36020k,0/ c c Machine constants for the Harris 220. c c data mcheps(1) / '20000000, '00000353 / c data minmag(1) / '20000000, '00000201 / c data maxmag(1) / '37777777, '00000177 / c c Machine constants for the Cray-1. c c data rmach(1) / 0377224000000000000000b / c data rmach(2) / 0200034000000000000000b / c data rmach(3) / 0577777777777777777776b / c c Machine constants for the Prime 400. c c data mcheps(1) / :10000000153 / c data minmag(1) / :10000000000 / c data maxmag(1) / :17777777777 / c c Machine constants for the VAX-11. c c data mcheps(1) / 13568 / c data minmag(1) / 128 / c data maxmag(1) / -32769 / c c Machine constants for IEEE machines. c data rmach(1) /1.192091E-07/ data rmach(2) /1.175495E-38/ data rmach(3) /3.402822E+38/ c spmpar = rmach(i) return c c Last card of function spmpar. c end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc REAL FUNCTION sgamma(a) C**********************************************************************C C C C C C (STANDARD-) G A M M A DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C PARAMETER A >= 1.0 ! C C C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C GENERATING GAMMA VARIATES BY A C C MODIFIED REJECTION TECHNIQUE. C C COMM. ACM, 25,1 (JAN. 1982), 47 - 54. C C C C STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER C C (STRAIGHTFORWARD IMPLEMENTATION) C C C C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C C SUNIF. The argument IR thus goes away. C C C C**********************************************************************C C C C PARAMETER 0.0 < A < 1.0 ! C C C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C COMPUTER METHODS FOR SAMPLING FROM GAMMA, C C BETA, POISSON AND BINOMIAL DISTRIBUTIONS. C C COMPUTING, 12 (1974), 223 - 246. C C C C (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER) C C C C**********************************************************************C C C C INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION C OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION C C COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K)) C COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K) C COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K) C DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121, + -.00007388,.00024511,.00024240/ DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921, + .1423657,-.1367177,.1233795/ DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/ C C PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A" C SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380 C DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/ C C SAVE STATEMENTS C SAVE aa,aaa,s2,s,d,q0,b,si,c C IF (a.EQ.aa) GO TO 10 IF (a.LT.1.0) GO TO 140 C C STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED C aa = a s2 = a - 0.5 s = sqrt(s2) d = sqrt32 - 12.0*s C C STEP 2: T=STANDARD NORMAL DEVIATE, C X=(S,1/2)-NORMAL DEVIATE. C IMMEDIATE ACCEPTANCE (I) C 10 t = snorm() x = s + 0.5*t sgamma = x*x IF (t.GE.0.0) RETURN C C STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S) C u = ranf() IF (d*u.LE.t*t*t) RETURN C C STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY C IF (a.EQ.aaa) GO TO 40 aaa = a r = 1.0/a q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r C C APPROXIMATION DEPENDING ON SIZE OF PARAMETER A C THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND C C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS C IF (a.LE.3.686) GO TO 30 IF (a.LE.13.022) GO TO 20 C C CASE 3: A .GT. 13.022 C b = 1.77 si = .75 c = .1515/s GO TO 40 C C CASE 2: 3.686 .LT. A .LE. 13.022 C 20 b = 1.654 + .0076*s2 si = 1.68/s + .275 c = .062/s + .024 GO TO 40 C C CASE 1: A .LE. 3.686 C 30 b = .463 + s + .178*s2 si = 1.235 c = .195/s - .079 + .16*s C C STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE C 40 IF (x.LE.0.0) GO TO 70 C C STEP 6: CALCULATION OF V AND QUOTIENT Q C v = t/ (s+s) IF (abs(v).LE.0.25) GO TO 50 q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v) GO TO 60 50 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v C C STEP 7: QUOTIENT ACCEPTANCE (Q) C 60 IF (alog(1.0-u).LE.q) RETURN C C STEP 8: E=STANDARD EXPONENTIAL DEVIATE C U= 0,1 -UNIFORM DEVIATE C T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE C 70 e = sexpo() u = ranf() u = u + u - 1.0 t = b + sign(si*e,u) IF (.NOT. (u.GE.0.0)) GO TO 80 t = b + si*e GO TO 90 80 t = b - si*e C C STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719 C 90 IF (t.LT. (-.7187449)) GO TO 70 C C STEP 10: CALCULATION OF V AND QUOTIENT Q C v = t/ (s+s) IF (abs(v).LE.0.25) GO TO 100 q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v) GO TO 110 100 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v C C STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8) C 110 IF (q.LE.0.0) GO TO 70 IF (q.LE.0.5) GO TO 120 w = exp(q) - 1.0 GO TO 130 120 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q C C IF T IS REJECTED, SAMPLE AGAIN AT STEP 8 C 130 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70 x = s + 0.5*t sgamma = x*x RETURN C C ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.)) C 140 aa = 0.0 b = 1.0 + .3678794*a 150 p = b*ranf() IF (p.GE.1.0) GO TO 160 sgamma = exp(alog(p)/a) IF (sexpo().LT.sgamma) GO TO 150 RETURN 160 sgamma = -alog((b-p)/a) IF (sexpo().LT. (1.0-a)*alog(sgamma)) GO TO 150 RETURN END REAL FUNCTION snorm() C**********************************************************************C C C C C C (STANDARD-) N O R M A L DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM C C SAMPLING FROM THE NORMAL DISTRIBUTION. C C MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. C C C C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' C C (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C C C C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C C SUNIF. The argument IR thus goes away. C C C C**********************************************************************C C DIMENSION a(32),d(31),t(31),h(31) C C THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND C H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE C DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991, + .2372021,.2776904,.3186394,.3601299,.4022501,.4450965, + .4887764,.5334097,.5791322,.6260990,.6744898,.7245144, + .7764218,.8305109,.8871466,.9467818,1.009990,1.077516, + 1.150349,1.229859,1.318011,1.417797,1.534121,1.675940, + 1.862732,2.153875/ DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243, + .1899108,.1812252,.1736014,.1668419,.1607967,.1553497, + .1504094,.1459026,.1417700,.1379632,.1344418,.1311722, + .1281260,.1252791,.1226109,.1201036,.1177417,.1155119, + .1134023,.1114027,.1095039/ DATA t/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2, + .7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1, + .1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1, + .2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1, + .4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1, + .9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980, + .5847031/ DATA h/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1, + .4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1, + .4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1, + .4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1, + .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1, + .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016, + .7010474/ C 10 u = ranf() s = 0.0 IF (u.GT.0.5) s = 1.0 u = u + u - s 20 u = 32.0*u i = int(u) IF (i.EQ.32) i = 31 IF (i.EQ.0) GO TO 100 C C START CENTER C 30 ustar = u - float(i) aa = a(i) 40 IF (ustar.LE.t(i)) GO TO 60 w = (ustar-t(i))*h(i) C C EXIT (BOTH CASES) C 50 y = aa + w snorm = y IF (s.EQ.1.0) snorm = -y RETURN C C CENTER CONTINUED C 60 u = ranf() w = u* (a(i+1)-aa) tt = (0.5*w+aa)*w GO TO 80 70 tt = u ustar = ranf() 80 IF (ustar.GT.tt) GO TO 50 90 u = ranf() IF (ustar.GE.u) GO TO 70 ustar = ranf() GO TO 40 C C START TAIL C 100 i = 6 aa = a(32) GO TO 120 110 aa = aa + d(i) i = i + 1 120 u = u + u IF (u.LT.1.0) GO TO 110 130 u = u - 1.0 140 w = u*d(i) tt = (0.5*w+aa)*w GO TO 160 150 tt = u 160 ustar = ranf() IF (ustar.GT.tt) GO TO 50 170 u = ranf() IF (ustar.GE.u) GO TO 150 u = ranf() GO TO 140 END REAL FUNCTION sexpo() C**********************************************************************C C C C C C (STANDARD-) E X P O N E N T I A L DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C COMPUTER METHODS FOR SAMPLING FROM THE C C EXPONENTIAL AND NORMAL DISTRIBUTIONS. C C COMM. ACM, 15,10 (OCT. 1972), 873 - 882. C C C C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM C C 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C C C C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C C SUNIF. The argument IR thus goes away. C C C C**********************************************************************C C DIMENSION q(8) EQUIVALENCE (q(1),q1) C C Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N C (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION C DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833, + .9999986,.9999999/ C 10 a = 0.0 u = ranf() GO TO 30 20 a = a + q1 30 u = u + u IF (u.LE.1.0) GO TO 20 40 u = u - 1.0 IF (u.GT.q1) GO TO 60 50 sexpo = a + u RETURN 60 i = 1 ustar = ranf() umin = ustar 70 ustar = ranf() IF (ustar.LT.umin) umin = ustar 80 i = i + 1 IF (u.GT.q(i)) GO TO 70 90 sexpo = a + umin*q1 RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc REAL FUNCTION ranf() C********************************************************************** C C REAL FUNCTION RANF() C RANDom number generator as a Function C C Returns a random floating point number from a uniform distribution C over 0 - 1 (endpoints of this interval are not returned) using the C current generator C C This is a transcription from Pascal to Fortran of routine C Uniform_01 from the paper C C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package C with Splitting Facilities." ACM Transactions on Mathematical C Software, 17:98-111 (1991) c c This was taken from the RANLIB library prepared by c Barry W. Brown and James Lovato c C C********************************************************************** C The functions ranf() and ignlgi() have been combined into ranf(). C********************************************************************** C C INTEGER FUNCTION IGNLGI() C GeNerate LarGe Integer C C Returns a random integer following a uniform distribution over C (1, 2147483562) using the current generator. C C This is a transcription from Pascal to Fortran of routine C Random from the paper C C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package C with Splitting Facilities." ACM Transactions on Mathematical C Software, 17:98-111 (1991) C C********************************************************************** C .. Parameters .. INTEGER numg PARAMETER (numg=32) C .. C .. Scalars in Common .. INTEGER sss1,sss2,a1,a1vw,a1w,a2,a2vw,a2w,m1,m2 C .. C .. Arrays in Common .. INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg), + lg2(numg) LOGICAL qanti(numg) C .. C .. Local Scalars .. INTEGER k,z C .. C .. C .. Common blocks .. COMMON /globe/m1,m2,sss1,sss2,a1,a2,a1w,a2w,a1vw,a2vw, + ig1,ig2,lg1,lg2,cg1, + cg2,qanti C .. C .. Save statement .. SAVE /globe/ C .. C .. Executable Statements .. C k = sss1/53668 sss1 = a1* (sss1-k*53668) - k*12211 IF (sss1.LT.0) sss1 = sss1 + m1 k = sss2/52774 sss2 = a2* (sss2-k*52774) - k*3791 IF (sss2.LT.0) sss2 = sss2 + m2 z = sss1 - sss2 IF (z.LT.1) z = z + m1 - 1 C C 4.656613057E-10 is 1/M1 M1 is set in a data statement in IGNLGI C and is currently 2147483563. If M1 changes, change this also. C ranf = z*4.656613057E-10 RETURN END LOGICAL FUNCTION qrgnin() C********************************************************************** C C LOGICAL FUNCTION QRGNIN() C Q Random GeNerators INitialized? C C A trivial routine to determine whether or not the random C number generator has been initialized. Returns .TRUE. if C it has, else .FALSE. C C********************************************************************** C .. Scalar Arguments .. LOGICAL qvalue C .. C .. Local Scalars .. LOGICAL qinit C .. C .. Entry Points .. LOGICAL qrgnsn C .. C .. Save statement .. SAVE qinit C .. C .. Data statements .. DATA qinit/.FALSE./ C .. C .. Executable Statements .. qrgnin = qinit RETURN ENTRY qrgnsn(qvalue) C********************************************************************** C C LOGICAL FUNCTION QRGNSN( QVALUE ) C Q Random GeNerators Set whether iNitialized C C Sets state of whether random number generator is initialized C to QVALUE. C C This routine is actually an entry in QRGNIN, hence it is a C logical function. It returns the (meaningless) value .TRUE. C C********************************************************************** qinit = qvalue qrgnsn = .TRUE. RETURN END SUBROUTINE getcgn(g) INTEGER g C********************************************************************** C C SUBROUTINE GETCGN(G) C Get GeNerator C C Returns in G the number of the current random number generator C C C Arguments C C C G <-- Number of the current random number generator (1..32) C INTEGER G C C********************************************************************** C INTEGER curntg,numg SAVE curntg PARAMETER (numg=32) DATA curntg/1/ C C .. C .. Scalars in Common .. INTEGER sss1,sss2,a1,a1vw,a1w,a2,a2vw,a2w,m1,m2 C .. C .. Arrays in Common .. INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg), + lg2(numg) LOGICAL qanti(numg) C .. C .. Common blocks .. COMMON /globe/m1,m2,sss1,sss2,a1,a2,a1w,a2w,a1vw,a2vw, + ig1,ig2,lg1,lg2,cg1, + cg2,qanti C .. C .. Save statement .. SAVE /globe/ C .. C .. Executable Statments g = curntg RETURN ENTRY setcgn(g) C********************************************************************** C C SUBROUTINE SETCGN( G ) C Set GeNerator C C Sets the current generator to G. All references to a generat C are to the current generator. C C C Arguments C C C G --> Number of the current random number generator (1..32) C INTEGER G C C********************************************************************** C C Abort if generator number out of range C 10 continue if((sss1.eq.-999).and.(sss2.eq.-999))then curntg = g else cg1(curntg)=sss1 cg2(curntg)=sss2 curntg = g sss1=cg1(curntg) sss2=cg2(curntg) endif RETURN END SUBROUTINE inrgcm() C********************************************************************** C C SUBROUTINE INRGCM() C INitialize Random number Generator CoMmon C C C Function C C C Initializes common area for random number generator. This saves C the nuisance of a BLOCK DATA routine and the difficulty of C assuring that the routine is loaded with the other routines. C C********************************************************************** C .. Parameters .. INTEGER numg PARAMETER (numg=32) C .. C .. Scalars in Common .. INTEGER sss1,sss2,a1,a1vw,a1w,a2,a2vw,a2w,m1,m2 C .. C .. Arrays in Common .. INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg), + lg2(numg) LOGICAL qanti(numg) C .. C .. Local Scalars .. INTEGER i LOGICAL qdum C .. C .. External Functions .. LOGICAL qrgnsn EXTERNAL qrgnsn C .. C .. Common blocks .. COMMON /globe/m1,m2,sss1,sss2,a1,a2,a1w,a2w,a1vw,a2vw, + ig1,ig2,lg1,lg2,cg1, + cg2,qanti C .. C .. Save statement .. SAVE /globe/ C .. C .. Executable Statements .. C V=20; W=30; C C A1W = MOD(A1**(2**W),M1) A2W = MOD(A2**(2**W),M2) C A1VW = MOD(A1**(2**(V+W)),M1) A2VW = MOD(A2**(2**(V+W)),M2) C C If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed. C An efficient way to precompute a**(2*j) MOD m is to start with C a and square it j times modulo m using the function MLTMOD. C m1 = 2147483563 m2 = 2147483399 a1 = 40014 a2 = 40692 a1w = 1033780774 a2w = 1494757890 a1vw = 2082007225 a2vw = 784306273 DO 10,i = 1,numg qanti(i) = .FALSE. 10 CONTINUE C C Tell the world that common has been initialized C qdum = qrgnsn(.TRUE.) RETURN END SUBROUTINE setall(iseed1,iseed2) C********************************************************************** C C SUBROUTINE SETALL(ISEED1,ISEED2) C SET ALL random number generators C C Sets the initial seed of generator 1 to ISEED1 and ISEED2. The C initial seeds of the other generators are set accordingly, and C all generators states are set to these seeds. C C This is a transcription from Pascal to Fortran of routine C Set_Initial_Seed from the paper C C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package C with Splitting Facilities." ACM Transactions on Mathematical C Software, 17:98-111 (1991) C C C Arguments C C C ISEED1 -> First of two integer seeds C INTEGER ISEED1 C C ISEED2 -> Second of two integer seeds C INTEGER ISEED1 C C********************************************************************** C .. Parameters .. INTEGER numg PARAMETER (numg=32) C .. C .. Scalar Arguments .. INTEGER iseed1,iseed2 LOGICAL qssd C .. C .. Scalars in Common .. INTEGER sss1,sss2,a1,a1vw,a1w,a2,a2vw,a2w,m1,m2 C .. C .. Arrays in Common .. INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg), + lg2(numg) LOGICAL qanti(numg) C .. C .. Local Scalars .. INTEGER g,ocgn LOGICAL qqssd C .. C .. External Functions .. INTEGER mltmod LOGICAL qrgnin EXTERNAL mltmod,qrgnin C .. C .. External Subroutines .. EXTERNAL getcgn,initgn,inrgcm,setcgn C .. C .. Common blocks .. COMMON /globe/m1,m2,sss1,sss2,a1,a2,a1w,a2w,a1vw,a2vw, $ ig1,ig2,lg1,lg2,cg1, + cg2,qanti C .. C .. Save statement .. SAVE /globe/,qqssd C .. C .. Data statements .. DATA qqssd/.FALSE./ C .. C .. Executable Statements .. sss1=-999 sss2=-999 C C TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE C HAS BEEN CALLED. C qqssd = .TRUE. CALL getcgn(ocgn) C C Initialize Common Block if Necessary C IF (.NOT. (qrgnin())) CALL inrgcm() ig1(1) = iseed1 ig2(1) = iseed2 CALL initgn(-1) DO 10,g = 2,numg ig1(g) = mltmod(a1vw,ig1(g-1),m1) ig2(g) = mltmod(a2vw,ig2(g-1),m2) CALL setcgn(g) CALL initgn(-1) 10 CONTINUE CALL setcgn(ocgn) sss1=cg1(ocgn) sss2=cg2(ocgn) RETURN ENTRY rgnqsd(qssd) C********************************************************************** C C SUBROUTINE RGNQSD C Random Number Generator Query SeeD set? C C Returns (LOGICAL) QSSD as .TRUE. if SETALL has been invoked, C otherwise returns .FALSE. C C********************************************************************** qssd = qqssd RETURN END INTEGER FUNCTION mltmod(a,s,m) C********************************************************************** C C INTEGER FUNCTION MLTMOD(A,S,M) C C Returns (A*S) MOD M C C This is a transcription from Pascal to Fortran of routine C MULtMod_Decompos from the paper C C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package C with Splitting Facilities." ACM Transactions on Mathematical C Software, 17:98-111 (1991) C C C Arguments C C C A, S, M --> C INTEGER A,S,M C C********************************************************************** C .. Parameters .. INTEGER h PARAMETER (h=32768) C .. C .. Scalar Arguments .. INTEGER a,m,s C .. C .. Local Scalars .. INTEGER a0,a1,k,p,q,qh,rh C .. C .. Executable Statements .. C C H = 2**((b-2)/2) where b = 32 because we are using a 32 bit C machine. On a different machine recompute H C 10 IF (.NOT. (a.LT.h)) GO TO 20 a0 = a p = 0 GO TO 120 20 a1 = a/h a0 = a - h*a1 qh = m/h rh = m - h*qh IF (.NOT. (a1.GE.h)) GO TO 50 a1 = a1 - h k = s/qh p = h* (s-k*qh) - k*rh 30 IF (.NOT. (p.LT.0)) GO TO 40 p = p + m GO TO 30 40 GO TO 60 50 p = 0 C C P = (A2*S*H)MOD M C 60 IF (.NOT. (a1.NE.0)) GO TO 90 q = m/a1 k = s/q p = p - k* (m-a1*q) IF (p.GT.0) p = p - m p = p + a1* (s-k*q) 70 IF (.NOT. (p.LT.0)) GO TO 80 p = p + m GO TO 70 80 CONTINUE 90 k = p/qh C C P = ((A2*H + A1)*S)MOD M C p = h* (p-k*qh) - k*rh 100 IF (.NOT. (p.LT.0)) GO TO 110 p = p + m GO TO 100 110 CONTINUE 120 IF (.NOT. (a0.NE.0)) GO TO 150 C C P = ((A2*H + A1)*H*S)MOD M C q = m/a0 k = s/q p = p - k* (m-a0*q) IF (p.GT.0) p = p - m p = p + a0* (s-k*q) 130 IF (.NOT. (p.LT.0)) GO TO 140 p = p + m GO TO 130 140 CONTINUE 150 mltmod = p C RETURN END SUBROUTINE initgn(isdtyp) C********************************************************************** C C SUBROUTINE INITGN(ISDTYP) C INIT-ialize current G-e-N-erator C C Reinitializes the state of the current generator C C This is a transcription from Pascal to Fortran of routine C Init_Generator from the paper C C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package C with Splitting Facilities." ACM Transactions on Mathematical C Software, 17:98-111 (1991) C C C Arguments C C C ISDTYP -> The state to which the generator is to be set C C ISDTYP = -1 => sets the seeds to their initial value C ISDTYP = 0 => sets the seeds to the first value of C the current block C ISDTYP = 1 => sets the seeds to the first value of C the next block C C INTEGER ISDTYP C C********************************************************************** C .. Parameters .. INTEGER numg PARAMETER (numg=32) C .. C .. Scalar Arguments .. INTEGER isdtyp C .. C .. Scalars in Common .. INTEGER sss1,sss2,a1,a1vw,a1w,a2,a2vw,a2w,m1,m2 C .. C .. Arrays in Common .. INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg), + lg2(numg) LOGICAL qanti(numg) C .. C .. Local Scalars .. INTEGER g C .. C .. External Functions .. LOGICAL qrgnin INTEGER mltmod EXTERNAL qrgnin,mltmod C .. C .. External Subroutines .. EXTERNAL getcgn C .. C .. Common blocks .. COMMON /globe/m1,m2,sss1,sss2,a1,a2,a1w,a2w,a1vw,a2vw, $ ig1,ig2,lg1,lg2,cg1, + cg2,qanti C .. C .. Save statement .. SAVE /globe/ C .. C .. Executable Statements .. C Abort unless random number generator initialized IF (qrgnin()) GO TO 10 10 CALL getcgn(g) IF ((-1).NE. (isdtyp)) GO TO 20 lg1(g) = ig1(g) lg2(g) = ig2(g) GO TO 50 20 IF ((0).NE. (isdtyp)) GO TO 30 CONTINUE GO TO 50 C do nothing 30 continue c 30 IF ((1).NE. (isdtyp)) GO TO 40 lg1(g) = mltmod(a1w,lg1(g),m1) lg2(g) = mltmod(a2w,lg2(g),m2) GO TO 50 50 cg1(g) = lg1(g) cg2(g) = lg2(g) RETURN END cccccccccccccccc stuff added to handle beta distributions ccccccc c written by Alfred H. Morris. C ALGORITHM 708, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 360-373z. C PROGRAM BTST (OUTPUT, TAPE6=OUTPUT) C----------------------------------------------------------------------- SUBROUTINE BRATIO (A, B, X, Y, W, W1, IERR) C----------------------------------------------------------------------- C C EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B) C C -------------------- C C IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1 C AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES C C W = IX(A,B) C W1 = 1 - IX(A,B) C C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. C IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND C W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED, C THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO C ONE OF THE FOLLOWING VALUES ... C C IERR = 1 IF A OR B IS NEGATIVE C IERR = 2 IF A = B = 0 C IERR = 3 IF X .LT. 0 OR X .GT. 1 C IERR = 4 IF Y .LT. 0 OR Y .GT. 1 C IERR = 5 IF X + Y .NE. 1 C IERR = 6 IF X = A = 0 C IERR = 7 IF Y = B = 0 C C-------------------- C WRITTEN BY ALFRED H. MORRIS, JR. C NAVAL SURFACE WARFARE CENTER C DAHLGREN, VIRGINIA C REVISED ... NOV 1991 C----------------------------------------------------------------------- REAL LAMBDA C----------------------------------------------------------------------- C C ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST C FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0 C EPS = SPMPAR(1) C C----------------------------------------------------------------------- W = 0.0 W1 = 0.0 IF (A .LT. 0.0 .OR. B .LT. 0.0) GO TO 300 IF (A .EQ. 0.0 .AND. B .EQ. 0.0) GO TO 310 IF (X .LT. 0.0 .OR. X .GT. 1.0) GO TO 320 IF (Y .LT. 0.0 .OR. Y .GT. 1.0) GO TO 330 Z = ((X + Y) - 0.5) - 0.5 IF (ABS(Z) .GT. 3.0*EPS) GO TO 340 C IERR = 0 IF (X .EQ. 0.0) GO TO 200 IF (Y .EQ. 0.0) GO TO 210 IF (A .EQ. 0.0) GO TO 211 IF (B .EQ. 0.0) GO TO 201 C EPS = AMAX1(EPS, 1.E-15) IF (AMAX1(A,B) .LT. 1.E-3*EPS) GO TO 230 C IND = 0 A0 = A B0 = B X0 = X Y0 = Y IF (AMIN1(A0, B0) .GT. 1.0) GO TO 30 C C PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1 C IF (X .LE. 0.5) GO TO 10 IND = 1 A0 = B B0 = A X0 = Y Y0 = X C 10 IF (B0 .LT. AMIN1(EPS,EPS*A0)) GO TO 80 IF (A0 .LT. AMIN1(EPS,EPS*B0) .AND. B0*X0 .LE. 1.0) GO TO 90 IF (AMAX1(A0, B0) .GT. 1.0) GO TO 20 IF (A0 .GE. AMIN1(0.2, B0)) GO TO 100 c IF (X0**A0 .LE. 0.9) GO TO 100 IF (exp(a0*log(x0)) .LE. 0.9) GO TO 100 IF (X0 .GE. 0.3) GO TO 110 N = 20 GO TO 130 C 20 IF (B0 .LE. 1.0) GO TO 100 IF (X0 .GE. 0.3) GO TO 110 IF (X0 .GE. 0.1) GO TO 21 c IF ((X0*B0)**A0 .LE. 0.7) GO TO 100 IF (exp(a0*log(x0*b0)) .LE. 0.7) GO TO 100 21 IF (B0 .GT. 15.0) GO TO 131 N = 20 GO TO 130 C C PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1 C 30 IF (A .GT. B) GO TO 31 LAMBDA = A - (A + B)*X GO TO 32 31 LAMBDA = (A + B)*Y - B 32 IF (LAMBDA .GE. 0.0) GO TO 40 IND = 1 A0 = B B0 = A X0 = Y Y0 = X LAMBDA = ABS(LAMBDA) C 40 IF (B0 .LT. 40.0 .AND. B0*X0 .LE. 0.7) GO TO 100 IF (B0 .LT. 40.0) GO TO 140 IF (A0 .GT. B0) GO TO 50 IF (A0 .LE. 100.0) GO TO 120 IF (LAMBDA .GT. 0.03*A0) GO TO 120 GO TO 180 50 IF (B0 .LE. 100.0) GO TO 120 IF (LAMBDA .GT. 0.03*B0) GO TO 120 GO TO 180 C C EVALUATION OF THE APPROPRIATE ALGORITHM C 80 W = FPSER(A0, B0, X0, EPS) W1 = 0.5 + (0.5 - W) GO TO 220 C 90 W1 = APSER(A0, B0, X0, EPS) W = 0.5 + (0.5 - W1) GO TO 220 C 100 W = BPSER(A0, B0, X0, EPS) W1 = 0.5 + (0.5 - W) GO TO 220 C 110 W1 = BPSER(B0, A0, Y0, EPS) W = 0.5 + (0.5 - W1) GO TO 220 C 120 W = BFRAC(A0, B0, X0, Y0, LAMBDA, 15.0*EPS) W1 = 0.5 + (0.5 - W) GO TO 220 C 130 W1 = BUP(B0, A0, Y0, X0, N, EPS) B0 = B0 + N 131 CALL BGRAT(B0, A0, Y0, X0, W1, 15.0*EPS, IERR1) W = 0.5 + (0.5 - W1) GO TO 220 C 140 N = B0 B0 = B0 - N IF (B0 .NE. 0.0) GO TO 141 N = N - 1 B0 = 1.0 141 W = BUP(B0, A0, Y0, X0, N, EPS) IF (X0 .GT. 0.7) GO TO 150 W = W + BPSER(A0, B0, X0, EPS) W1 = 0.5 + (0.5 - W) GO TO 220 C 150 IF (A0 .GT. 15.0) GO TO 151 N = 20 W = W + BUP(A0, B0, X0, Y0, N, EPS) A0 = A0 + N 151 CALL BGRAT(A0, B0, X0, Y0, W, 15.0*EPS, IERR1) W1 = 0.5 + (0.5 - W) GO TO 220 C 180 W = BASYM(A0, B0, LAMBDA, 100.0*EPS) W1 = 0.5 + (0.5 - W) GO TO 220 C C TERMINATION OF THE PROCEDURE C 200 IF (A .EQ. 0.0) GO TO 350 201 W = 0.0 W1 = 1.0 RETURN C 210 IF (B .EQ. 0.0) GO TO 360 211 W = 1.0 W1 = 0.0 RETURN C 220 IF (IND .EQ. 0) RETURN T = W W = W1 W1 = T RETURN C C PROCEDURE FOR A AND B .LT. 1.E-3*EPS C 230 W = B/(A + B) W1 = A/(A + B) RETURN C C ERROR RETURN C 300 IERR = 1 RETURN 310 IERR = 2 RETURN 320 IERR = 3 RETURN 330 IERR = 4 RETURN 340 IERR = 5 RETURN 350 IERR = 6 RETURN 360 IERR = 7 RETURN END REAL FUNCTION FPSER (A, B, X, EPS) C----------------------------------------------------------------------- C C EVALUATION OF I (A,B) C X C C FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5. C C----------------------------------------------------------------------- C C SET FPSER = X**A C FPSER = 1.0 IF (A .LE. 1.E-3*EPS) GO TO 10 FPSER = 0.0 T = A*ALOG(X) IF (T .LT. EXPARG(1)) RETURN FPSER = EXP(T) C C NOTE THAT 1/B(A,B) = B C 10 FPSER = (B/A)*FPSER TOL = EPS/A AN = A + 1.0 T = X S = T/AN 20 AN = AN + 1.0 T = X*T C = T/AN S = S + C IF (ABS(C) .GT. TOL) GO TO 20 C FPSER = FPSER*(1.0 + A*S) RETURN END REAL FUNCTION APSER (A, B, X, EPS) C----------------------------------------------------------------------- C APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR C A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN C A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED. C----------------------------------------------------------------------- REAL J C-------------------- DATA G/.577215664901533/ C-------------------- BX = B*X T = X - BX IF (B*EPS .GT. 2.E-2) GO TO 10 C = ALOG(X) + PSI(B) + G + T GO TO 20 10 C = ALOG(BX) + G + T C 20 TOL = 5.0*EPS*ABS(C) J = 1.0 S = 0.0 30 J = J + 1.0 T = T*(X - BX/J) AJ = T/J S = S + AJ IF (ABS(AJ) .GT. TOL) GO TO 30 C APSER = -A*(C + S) RETURN END REAL FUNCTION BPSER(A, B, X, EPS) C----------------------------------------------------------------------- C POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1 C OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED. C----------------------------------------------------------------------- REAL N C BPSER = 0.0 IF (X .EQ. 0.0) RETURN C----------------------------------------------------------------------- C COMPUTE THE FACTOR X**A/(A*BETA(A,B)) C----------------------------------------------------------------------- A0 = AMIN1(A,B) IF (A0 .LT. 1.0) GO TO 10 Z = A*ALOG(X) - BETALN(A,B) BPSER = EXP(Z)/A GO TO 70 10 B0 = AMAX1(A,B) IF (B0 .GE. 8.0) GO TO 60 IF (B0 .GT. 1.0) GO TO 40 C C PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1 C c Comment this out to get rid of pow_dd error when dyn.loading. c BPSER = X**A bpser=exp(a*log(x)) IF (BPSER .EQ. 0.0) RETURN C APB = A + B IF (APB .GT. 1.0) GO TO 20 Z = 1.0 + GAM1(APB) GO TO 30 20 U = DBLE(A) + DBLE(B) - 1.D0 Z = (1.0 + GAM1(U))/APB C 30 C = (1.0 + GAM1(A))*(1.0 + GAM1(B))/Z BPSER = BPSER*C*(B/APB) GO TO 70 C C PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8 C 40 U = GAMLN1(A0) M = B0 - 1.0 IF (M .LT. 1) GO TO 50 C = 1.0 DO 41 I = 1,M B0 = B0 - 1.0 41 C = C*(B0/(A0 + B0)) U = ALOG(C) + U C 50 Z = A*ALOG(X) - U B0 = B0 - 1.0 APB = A0 + B0 IF (APB .GT. 1.0) GO TO 51 T = 1.0 + GAM1(APB) GO TO 52 51 U = DBLE(A0) + DBLE(B0) - 1.D0 T = (1.0 + GAM1(U))/APB 52 BPSER = EXP(Z)*(A0/A)*(1.0 + GAM1(B0))/T GO TO 70 C C PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8 C 60 U = GAMLN1(A0) + ALGDIV(A0,B0) Z = A*ALOG(X) - U BPSER = (A0/A)*EXP(Z) 70 IF (BPSER .EQ. 0.0 .OR. A .LE. 0.1*EPS) RETURN C----------------------------------------------------------------------- C COMPUTE THE SERIES C----------------------------------------------------------------------- SUM = 0.0 N = 0.0 C = 1.0 TOL = EPS/A 100 N = N + 1.0 C = C*(0.5 + (0.5 - B/N))*X W = C/(A + N) SUM = SUM + W IF (ABS(W) .GT. TOL) GO TO 100 BPSER = BPSER*(1.0 + A*SUM) RETURN END REAL FUNCTION BUP(A, B, X, Y, N, EPS) C----------------------------------------------------------------------- C EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER. C EPS IS THE TOLERANCE USED. C----------------------------------------------------------------------- REAL L C C OBTAIN THE SCALING FACTOR EXP(-MU) AND C EXP(MU)*(X**A*Y**B/BETA(A,B))/A C APB = A + B AP1 = A + 1.0 MU = 0 D = 1.0 IF (N .EQ. 1 .OR. A .LT. 1.0) GO TO 10 IF (APB .LT. 1.1*AP1) GO TO 10 MU = ABS(EXPARG(1)) K = EXPARG(0) IF (K .LT. MU) MU = K T = MU D = EXP(-T) C 10 BUP = BRCMP1(MU,A,B,X,Y)/A IF (N .EQ. 1 .OR. BUP .EQ. 0.0) RETURN NM1 = N - 1 W = D C C LET K BE THE INDEX OF THE MAXIMUM TERM C K = 0 IF (B .LE. 1.0) GO TO 40 IF (Y .GT. 1.E-4) GO TO 20 K = NM1 GO TO 30 20 R = (B - 1.0)*X/Y - A IF (R .LT. 1.0) GO TO 40 K = NM1 T = NM1 IF (R .LT. T) K = R C C ADD THE INCREASING TERMS OF THE SERIES C 30 DO 31 I = 1,K L = I - 1 D = ((APB + L)/(AP1 + L))*X*D W = W + D 31 CONTINUE IF (K .EQ. NM1) GO TO 50 C C ADD THE REMAINING TERMS OF THE SERIES C 40 KP1 = K + 1 DO 41 I = KP1,NM1 L = I - 1 D = ((APB + L)/(AP1 + L))*X*D W = W + D IF (D .LE. EPS*W) GO TO 50 41 CONTINUE C C TERMINATE THE PROCEDURE C 50 BUP = BUP*W RETURN END REAL FUNCTION BFRAC(A, B, X, Y, LAMBDA, EPS) C----------------------------------------------------------------------- C CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1. C IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B. C----------------------------------------------------------------------- REAL LAMBDA, N C-------------------- BFRAC = BRCOMP(A,B,X,Y) IF (BFRAC .EQ. 0.0) RETURN C C = 1.0 + LAMBDA C0 = B/A C1 = 1.0 + 1.0/A YP1 = Y + 1.0 C N = 0.0 P = 1.0 S = A + 1.0 AN = 0.0 BN = 1.0 ANP1 = 1.0 BNP1 = C/C1 R = C1/C C C CONTINUED FRACTION CALCULATION C 10 N = N + 1.0 T = N/A W = N*(B - N)*X E = A/S ALPHA = (P*(P + C0)*E*E)*(W*X) E = (1.0 + T)/(C1 + T + T) BETA = N + W/S + E*(C + N*YP1) P = 1.0 + T S = S + 2.0 C C UPDATE AN, BN, ANP1, AND BNP1 C T = ALPHA*AN + BETA*ANP1 AN = ANP1 ANP1 = T T = ALPHA*BN + BETA*BNP1 BN = BNP1 BNP1 = T C R0 = R R = ANP1/BNP1 IF (ABS(R - R0) .LE. EPS*R) GO TO 20 C C RESCALE AN, BN, ANP1, AND BNP1 C AN = AN/BNP1 BN = BN/BNP1 ANP1 = R BNP1 = 1.0 GO TO 10 C C TERMINATION C 20 BFRAC = BFRAC*R RETURN END REAL FUNCTION BRCOMP (A, B, X, Y) C----------------------------------------------------------------------- C EVALUATION OF X**A*Y**B/BETA(A,B) C----------------------------------------------------------------------- REAL LAMBDA, LNX, LNY C----------------- C CONST = 1/SQRT(2*PI) C----------------- DATA CONST/.398942280401433/ C BRCOMP = 0.0 IF (X .EQ. 0.0 .OR. Y .EQ. 0.0) RETURN A0 = AMIN1(A,B) IF (A0 .GE. 8.0) GO TO 100 C IF (X .GT. 0.375) GO TO 10 LNX = ALOG(X) LNY = ALNREL(-X) GO TO 20 10 IF (Y .GT. 0.375) GO TO 11 LNX = ALNREL(-Y) LNY = ALOG(Y) GO TO 20 11 LNX = ALOG(X) LNY = ALOG(Y) C 20 Z = A*LNX + B*LNY IF (A0 .LT. 1.0) GO TO 30 Z = Z - BETALN(A,B) BRCOMP = EXP(Z) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .LT. 1 OR B .LT. 1 C----------------------------------------------------------------------- 30 B0 = AMAX1(A,B) IF (B0 .GE. 8.0) GO TO 80 IF (B0 .GT. 1.0) GO TO 60 C C ALGORITHM FOR B0 .LE. 1 C BRCOMP = EXP(Z) IF (BRCOMP .EQ. 0.0) RETURN C APB = A + B IF (APB .GT. 1.0) GO TO 40 Z = 1.0 + GAM1(APB) GO TO 50 40 U = DBLE(A) + DBLE(B) - 1.D0 Z = (1.0 + GAM1(U))/APB C 50 C = (1.0 + GAM1(A))*(1.0 + GAM1(B))/Z BRCOMP = BRCOMP*(A0*C)/(1.0 + A0/B0) RETURN C C ALGORITHM FOR 1 .LT. B0 .LT. 8 C 60 U = GAMLN1(A0) N = B0 - 1.0 IF (N .LT. 1) GO TO 70 C = 1.0 DO 61 I = 1,N B0 = B0 - 1.0 C = C*(B0/(A0 + B0)) 61 CONTINUE U = ALOG(C) + U C 70 Z = Z - U B0 = B0 - 1.0 APB = A0 + B0 IF (APB .GT. 1.0) GO TO 71 T = 1.0 + GAM1(APB) GO TO 72 71 U = DBLE(A0) + DBLE(B0) - 1.D0 T = (1.0 + GAM1(U))/APB 72 BRCOMP = A0*EXP(Z)*(1.0 + GAM1(B0))/T RETURN C C ALGORITHM FOR B0 .GE. 8 C 80 U = GAMLN1(A0) + ALGDIV(A0,B0) BRCOMP = A0*EXP(Z - U) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .GE. 8 AND B .GE. 8 C----------------------------------------------------------------------- 100 IF (A .GT. B) GO TO 101 H = A/B X0 = H/(1.0 + H) Y0 = 1.0/(1.0 + H) LAMBDA = A - (A + B)*X GO TO 110 101 H = B/A X0 = 1.0/(1.0 + H) Y0 = H/(1.0 + H) LAMBDA = (A + B)*Y - B C 110 E = -LAMBDA/A IF (ABS(E) .GT. 0.6) GO TO 111 U = RLOG1(E) GO TO 120 111 U = E - ALOG(X/X0) C 120 E = LAMBDA/B IF (ABS(E) .GT. 0.6) GO TO 121 V = RLOG1(E) GO TO 130 121 V = E - ALOG(Y/Y0) C 130 Z = EXP(-(A*U + B*V)) BRCOMP = CONST*SQRT(B*X0)*Z*EXP(-BCORR(A,B)) RETURN END REAL FUNCTION BRCMP1 (MU, A, B, X, Y) C----------------------------------------------------------------------- C EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B)) C----------------------------------------------------------------------- REAL LAMBDA, LNX, LNY C----------------- C CONST = 1/SQRT(2*PI) C----------------- DATA CONST/.398942280401433/ C A0 = AMIN1(A,B) IF (A0 .GE. 8.0) GO TO 100 C IF (X .GT. 0.375) GO TO 10 LNX = ALOG(X) LNY = ALNREL(-X) GO TO 20 10 IF (Y .GT. 0.375) GO TO 11 LNX = ALNREL(-Y) LNY = ALOG(Y) GO TO 20 11 LNX = ALOG(X) LNY = ALOG(Y) C 20 Z = A*LNX + B*LNY IF (A0 .LT. 1.0) GO TO 30 Z = Z - BETALN(A,B) BRCMP1 = ESUM(MU,Z) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .LT. 1 OR B .LT. 1 C----------------------------------------------------------------------- 30 B0 = AMAX1(A,B) IF (B0 .GE. 8.0) GO TO 80 IF (B0 .GT. 1.0) GO TO 60 C C ALGORITHM FOR B0 .LE. 1 C BRCMP1 = ESUM(MU,Z) IF (BRCMP1 .EQ. 0.0) RETURN C APB = A + B IF (APB .GT. 1.0) GO TO 40 Z = 1.0 + GAM1(APB) GO TO 50 40 U = DBLE(A) + DBLE(B) - 1.D0 Z = (1.0 + GAM1(U))/APB C 50 C = (1.0 + GAM1(A))*(1.0 + GAM1(B))/Z BRCMP1 = BRCMP1*(A0*C)/(1.0 + A0/B0) RETURN C C ALGORITHM FOR 1 .LT. B0 .LT. 8 C 60 U = GAMLN1(A0) N = B0 - 1.0 IF (N .LT. 1) GO TO 70 C = 1.0 DO 61 I = 1,N B0 = B0 - 1.0 C = C*(B0/(A0 + B0)) 61 CONTINUE U = ALOG(C) + U C 70 Z = Z - U B0 = B0 - 1.0 APB = A0 + B0 IF (APB .GT. 1.0) GO TO 71 T = 1.0 + GAM1(APB) GO TO 72 71 U = DBLE(A0) + DBLE(B0) - 1.D0 T = (1.0 + GAM1(U))/APB 72 BRCMP1 = A0*ESUM(MU,Z)*(1.0 + GAM1(B0))/T RETURN C C ALGORITHM FOR B0 .GE. 8 C 80 U = GAMLN1(A0) + ALGDIV(A0,B0) BRCMP1 = A0*ESUM(MU,Z - U) RETURN C----------------------------------------------------------------------- C PROCEDURE FOR A .GE. 8 AND B .GE. 8 C----------------------------------------------------------------------- 100 IF (A .GT. B) GO TO 101 H = A/B X0 = H/(1.0 + H) Y0 = 1.0/(1.0 + H) LAMBDA = A - (A + B)*X GO TO 110 101 H = B/A X0 = 1.0/(1.0 + H) Y0 = H/(1.0 + H) LAMBDA = (A + B)*Y - B C 110 E = -LAMBDA/A IF (ABS(E) .GT. 0.6) GO TO 111 U = RLOG1(E) GO TO 120 111 U = E - ALOG(X/X0) C 120 E = LAMBDA/B IF (ABS(E) .GT. 0.6) GO TO 121 V = RLOG1(E) GO TO 130 121 V = E - ALOG(Y/Y0) C 130 Z = ESUM(MU,-(A*U + B*V)) BRCMP1 = CONST*SQRT(B*X0)*Z*EXP(-BCORR(A,B)) RETURN END SUBROUTINE BGRAT(A, B, X, Y, W, EPS, IERR) C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B. C THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED C THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED. C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS. C----------------------------------------------------------------------- REAL J, L, LNX, NU, N2 REAL C(30), D(30) C BM1 = (B - 0.5) - 0.5 NU = A + 0.5*BM1 IF (Y .GT. 0.375) GO TO 10 LNX = ALNREL(-Y) GO TO 11 10 LNX = ALOG(X) 11 Z = -NU*LNX IF (B*Z .EQ. 0.0) GO TO 100 C C COMPUTATION OF THE EXPANSION C SET R = EXP(-Z)*Z**B/GAMMA(B) C R = B*(1.0 + GAM1(B))*EXP(B*ALOG(Z)) R = R*EXP(A*LNX)*EXP(0.5*BM1*LNX) U = ALGDIV(B,A) + B*ALOG(NU) U = R*EXP(-U) IF (U .EQ. 0.0) GO TO 100 CALL GRAT1(B,Z,R,P,Q,EPS) C V = 0.25*(1.0/NU)**2 T2 = 0.25*LNX*LNX L = W/U J = Q/R SUM = J T = 1.0 CN = 1.0 N2 = 0.0 DO 22 N = 1,30 BP2N = B + N2 J = (BP2N*(BP2N + 1.0)*J + (Z + BP2N + 1.0)*T)*V N2 = N2 + 2.0 T = T*T2 CN = CN/(N2*(N2 + 1.0)) C(N) = CN S = 0.0 IF (N .EQ. 1) GO TO 21 NM1 = N - 1 COEF = B - N DO 20 I = 1,NM1 S = S + COEF*C(I)*D(N-I) 20 COEF = COEF + B 21 D(N) = BM1*CN + S/N DJ = D(N)*J SUM = SUM + DJ IF (SUM .LE. 0.0) GO TO 100 IF (ABS(DJ) .LE. EPS*(SUM + L)) GO TO 30 22 CONTINUE C C ADD THE RESULTS TO W C 30 IERR = 0 W = W + U*SUM RETURN C C THE EXPANSION CANNOT BE COMPUTED C 100 IERR = 1 RETURN END SUBROUTINE GRAT1 (A,X,R,P,Q,EPS) REAL J, L C----------------------------------------------------------------------- C EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS C P(A,X) AND Q(A,X) C C IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED. C THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A). C----------------------------------------------------------------------- IF (A*X .EQ. 0.0) GO TO 130 IF (A .EQ. 0.5) GO TO 120 IF (X .LT. 1.1) GO TO 10 GO TO 50 C C TAYLOR SERIES FOR P(A,X)/X**A C 10 AN = 3.0 C = X SUM = X/(A + 3.0) TOL = 0.1*EPS/(A + 1.0) 11 AN = AN + 1.0 C = -C*(X/AN) T = C/(A + AN) SUM = SUM + T IF (ABS(T) .GT. TOL) GO TO 11 J = A*X*((SUM/6.0 - 0.5/(A + 2.0))*X + 1.0/(A + 1.0)) C Z = A*ALOG(X) H = GAM1(A) G = 1.0 + H IF (X .LT. 0.25) GO TO 20 IF (A .LT. X/2.59) GO TO 40 GO TO 30 20 IF (Z .GT. -.13394) GO TO 40 C 30 W = EXP(Z) P = W*G*(0.5 + (0.5 - J)) Q = 0.5 + (0.5 - P) RETURN C 40 L = REXP(Z) W = 0.5 + (0.5 + L) Q = (W*J - L)*G - H IF (Q .LT. 0.0) GO TO 110 P = 0.5 + (0.5 - Q) RETURN C C CONTINUED FRACTION EXPANSION C 50 A2NM1 = 1.0 A2N = 1.0 B2NM1 = X B2N = X + (1.0 - A) C = 1.0 51 A2NM1 = X*A2N + C*A2NM1 B2NM1 = X*B2N + C*B2NM1 AM0 = A2NM1/B2NM1 C = C + 1.0 CMA = C - A A2N = A2NM1 + CMA*A2N B2N = B2NM1 + CMA*B2N AN0 = A2N/B2N IF (ABS(AN0 - AM0) .GE. EPS*AN0) GO TO 51 Q = R*AN0 P = 0.5 + (0.5 - Q) RETURN C C SPECIAL CASES C 100 P = 0.0 Q = 1.0 RETURN C 110 P = 1.0 Q = 0.0 RETURN C 120 IF (X .GE. 0.25) GO TO 121 P = ERFff(SQRT(X)) Q = 0.5 + (0.5 - P) RETURN 121 Q = ERFC1(0,SQRT(X)) P = 0.5 + (0.5 - Q) RETURN C 130 IF (X .LE. A) GO TO 100 GO TO 110 END REAL FUNCTION BASYM(A, B, LAMBDA, EPS) C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B. C LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED. C IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT C A AND B ARE GREATER THAN OR EQUAL TO 15. C----------------------------------------------------------------------- REAL J0, J1, LAMBDA REAL A0(21), B0(21), C(21), D(21) C------------------------ C ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP C ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN. C THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1. C DATA NUM/20/ C------------------------ C E0 = 2/SQRT(PI) C E1 = 2**(-3/2) C------------------------ DATA E0/1.12837916709551/, E1/.353553390593274/ C------------------------ BASYM = 0.0 IF (A .GE. B) GO TO 10 H = A/B R0 = 1.0/(1.0 + H) R1 = (B - A)/B W0 = 1.0/SQRT(A*(1.0 + H)) GO TO 20 10 H = B/A R0 = 1.0/(1.0 + H) R1 = (B - A)/A W0 = 1.0/SQRT(B*(1.0 + H)) C 20 F = A*RLOG1(-LAMBDA/A) + B*RLOG1(LAMBDA/B) T = EXP(-F) IF (T .EQ. 0.0) RETURN Z0 = SQRT(F) Z = 0.5*(Z0/E1) Z2 = F + F C A0(1) = (2.0/3.0)*R1 C(1) = - 0.5*A0(1) D(1) = - C(1) J0 = (0.5/E0)*ERFC1(1,Z0) J1 = E1 SUM = J0 + D(1)*W0*J1 C S = 1.0 H2 = H*H HN = 1.0 W = W0 ZNM1 = Z ZN = Z2 DO 50 N = 2, NUM, 2 HN = H2*HN A0(N) = 2.0*R0*(1.0 + H*HN)/(N + 2.0) NP1 = N + 1 S = S + HN A0(NP1) = 2.0*R1*S/(N + 3.0) C DO 41 I = N, NP1 R = -0.5*(I + 1.0) B0(1) = R*A0(1) DO 31 M = 2, I BSUM = 0.0 MM1 = M - 1 DO 30 J = 1, MM1 MMJ = M - J 30 BSUM = BSUM + (J*R - MMJ)*A0(J)*B0(MMJ) 31 B0(M) = R*A0(M) + BSUM/M C(I) = B0(I)/(I + 1.0) C DSUM = 0.0 IM1 = I - 1 DO 40 J = 1, IM1 IMJ = I - J 40 DSUM = DSUM + D(IMJ)*C(J) 41 D(I) = -(DSUM + C(I)) C J0 = E1*ZNM1 + (N - 1.0)*J0 J1 = E1*ZN + N*J1 ZNM1 = Z2*ZNM1 ZN = Z2*ZN W = W0*W T0 = D(N)*W*J0 W = W0*W T1 = D(NP1)*W*J1 SUM = SUM + (T0 + T1) IF ((ABS(T0) + ABS(T1)) .LE. EPS*SUM) GO TO 60 50 CONTINUE C 60 U = EXP(-BCORR(A,B)) BASYM = E0*T*U*SUM RETURN END REAL FUNCTION ESUM (MU, X) C----------------------------------------------------------------------- C EVALUATION OF EXP(MU + X) C----------------------------------------------------------------------- IF (X .GT. 0.0) GO TO 10 C IF (MU .LT. 0) GO TO 20 W = MU + X IF (W .GT. 0.0) GO TO 20 ESUM = EXP(W) RETURN C 10 IF (MU .GT. 0) GO TO 20 W = MU + X IF (W .LT. 0.0) GO TO 20 ESUM = EXP(W) RETURN C 20 W = MU ESUM = EXP(W)*EXP(X) RETURN END REAL FUNCTION RLOG1(X) C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION X - LN(1 + X) C----------------------------------------------------------------------- DATA A/.566749439387324E-01/ DATA B/.456512608815524E-01/ C------------------------ DATA P0/ .333333333333333E+00/, P1/-.224696413112536E+00/, * P2/ .620886815375787E-02/ DATA Q1/-.127408923933623E+01/, Q2/ .354508718369557E+00/ C------------------------ IF (X .LT. -0.39 .OR. X .GT. 0.57) GO TO 100 IF (X .LT. -0.18) GO TO 10 IF (X .GT. 0.18) GO TO 20 C C ARGUMENT REDUCTION C H = X W1 = 0.0 GO TO 30 C 10 H = DBLE(X) + 0.3D0 H = H/0.7 W1 = A - H*0.3 GO TO 30 C 20 H = 0.75D0*DBLE(X) - 0.25D0 W1 = B + H/3.0 C C SERIES EXPANSION C 30 R = H/(H + 2.0) T = R*R W = ((P2*T + P1)*T + P0)/((Q2*T + Q1)*T + 1.0) RLOG1 = 2.0*T*(1.0/(1.0 - R) - R*W) + W1 RETURN C C 100 W = (X + 0.5) + 0.5 RLOG1 = X - ALOG(W) RETURN END REAL FUNCTION PSI(XX) C--------------------------------------------------------------------- C C EVALUATION OF THE DIGAMMA FUNCTION C C ----------- C C PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT C BE COMPUTED. C C THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV C APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY C CODY, STRECOK AND THACHER. C C--------------------------------------------------------------------- C PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK C PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY C A.H. MORRIS (NSWC). C--------------------------------------------------------------------- c change made be me so I can get rid of the ipmpar function. parameter(xipm3=2147483647.0) REAL P1(7), P2(4), Q1(6), Q2(4) DOUBLE PRECISION DX0 C--------------------------------------------------------------------- C C PIOV4 = PI/4 C DX0 = ZERO OF PSI TO EXTENDED PRECISION C C--------------------------------------------------------------------- DATA PIOV4/.785398163397448E0/ DATA DX0/1.461632144968362341262659542325721325D0/ C--------------------------------------------------------------------- C C COEFFICIENTS FOR RATIONAL APPROXIMATION OF C PSI(X) / (X - X0), 0.5 .LE. X .LE. 3.0 C C--------------------------------------------------------------------- DATA P1(1)/.895385022981970E-02/, P1(2)/.477762828042627E+01/, * P1(3)/.142441585084029E+03/, P1(4)/.118645200713425E+04/, * P1(5)/.363351846806499E+04/, P1(6)/.413810161269013E+04/, * P1(7)/.130560269827897E+04/ DATA Q1(1)/.448452573429826E+02/, Q1(2)/.520752771467162E+03/, * Q1(3)/.221000799247830E+04/, Q1(4)/.364127349079381E+04/, * Q1(5)/.190831076596300E+04/, Q1(6)/.691091682714533E-05/ C--------------------------------------------------------------------- C C COEFFICIENTS FOR RATIONAL APPROXIMATION OF C PSI(X) - LN(X) + 1 / (2*X), X .GT. 3.0 C C--------------------------------------------------------------------- DATA P2(1)/-.212940445131011E+01/, P2(2)/-.701677227766759E+01/, * P2(3)/-.448616543918019E+01/, P2(4)/-.648157123766197E+00/ DATA Q2(1)/ .322703493791143E+02/, Q2(2)/ .892920700481861E+02/, * Q2(3)/ .546117738103215E+02/, Q2(4)/ .777788548522962E+01/ C--------------------------------------------------------------------- C C MACHINE DEPENDENT CONSTANTS ... C C XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT C WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED C AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE C ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH C PSI MAY BE REPRESENTED AS ALOG(X). C C XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X) C MAY BE REPRESENTED BY 1/X. C C--------------------------------------------------------------------- c commenting this out so I can get rid of ipmpar. c XMAX1 = IPMPAR(3) XMAX1 = AMIN1(xipm3, 1.0/SPMPAR(1)) XSMALL = 1.E-9 C--------------------------------------------------------------------- X = XX AUG = 0.0E0 IF (X .GE. 0.5E0) GO TO 200 C--------------------------------------------------------------------- C X .LT. 0.5, USE REFLECTION FORMULA C PSI(1-X) = PSI(X) + PI * COTAN(PI*X) C--------------------------------------------------------------------- IF (ABS(X) .GT. XSMALL) GO TO 100 IF (X .EQ. 0.0E0) GO TO 400 C--------------------------------------------------------------------- C 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE C FOR PI*COTAN(PI*X) C--------------------------------------------------------------------- AUG = -1.0E0 / X GO TO 150 C--------------------------------------------------------------------- C REDUCTION OF ARGUMENT FOR COTAN C--------------------------------------------------------------------- 100 W = - X SGN = PIOV4 IF (W .GT. 0.0E0) GO TO 120 W = - W SGN = -SGN C--------------------------------------------------------------------- C MAKE AN ERROR EXIT IF X .LE. -XMAX1 C--------------------------------------------------------------------- 120 IF (W .GE. XMAX1) GO TO 400 NQ = INT(W) W = W - FLOAT(NQ) NQ = INT(W*4.0E0) W = 4.0E0 * (W - FLOAT(NQ) * .25E0) C--------------------------------------------------------------------- C W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X. C ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST C QUADRANT AND DETERMINE SIGN C--------------------------------------------------------------------- N = NQ / 2 IF ((N+N) .NE. NQ) W = 1.0E0 - W Z = PIOV4 * W M = N / 2 IF ((M+M) .NE. N) SGN = - SGN C--------------------------------------------------------------------- C DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X) C--------------------------------------------------------------------- N = (NQ + 1) / 2 M = N / 2 M = M + M IF (M .NE. N) GO TO 140 C--------------------------------------------------------------------- C CHECK FOR SINGULARITY C--------------------------------------------------------------------- IF (Z .EQ. 0.0E0) GO TO 400 C--------------------------------------------------------------------- C USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND C SIN/COS AS A SUBSTITUTE FOR TAN C--------------------------------------------------------------------- AUG = SGN * ((COS(Z) / SIN(Z)) * 4.0E0) GO TO 150 140 AUG = SGN * ((SIN(Z) / COS(Z)) * 4.0E0) 150 X = 1.0E0 - X 200 IF (X .GT. 3.0E0) GO TO 300 C--------------------------------------------------------------------- C 0.5 .LE. X .LE. 3.0 C--------------------------------------------------------------------- DEN = X UPPER = P1(1) * X C DO 210 I = 1, 5 DEN = (DEN + Q1(I)) * X UPPER = (UPPER + P1(I+1)) * X 210 CONTINUE C DEN = (UPPER + P1(7)) / (DEN + Q1(6)) XMX0 = DBLE(X) - DX0 PSI = DEN * XMX0 + AUG RETURN C--------------------------------------------------------------------- C IF X .GE. XMAX1, PSI = LN(X) C--------------------------------------------------------------------- 300 IF (X .GE. XMAX1) GO TO 350 C--------------------------------------------------------------------- C 3.0 .LT. X .LT. XMAX1 C--------------------------------------------------------------------- W = 1.0E0 / (X * X) DEN = W UPPER = P2(1) * W C DO 310 I = 1, 3 DEN = (DEN + Q2(I)) * W UPPER = (UPPER + P2(I+1)) * W 310 CONTINUE C AUG = UPPER / (DEN + Q2(4)) - 0.5E0 / X + AUG 350 PSI = AUG + ALOG(X) RETURN C--------------------------------------------------------------------- C ERROR RETURN C--------------------------------------------------------------------- 400 PSI = 0.0E0 RETURN END REAL FUNCTION BETALN (A0, B0) C----------------------------------------------------------------------- C EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION C----------------------------------------------------------------------- C E = 0.5*LN(2*PI) C-------------------------- DATA E /.918938533204673/ C-------------------------- A = AMIN1(A0,B0) B = AMAX1(A0,B0) IF (A .GE. 8.0) GO TO 60 IF (A .GE. 1.0) GO TO 20 C----------------------------------------------------------------------- C PROCEDURE WHEN A .LT. 1 C----------------------------------------------------------------------- IF (B .GE. 8.0) GO TO 10 BETALN = GAMLN(A) + (GAMLN(B) - GAMLN(A + B)) RETURN 10 BETALN = GAMLN(A) + ALGDIV(A,B) RETURN C----------------------------------------------------------------------- C PROCEDURE WHEN 1 .LE. A .LT. 8 C----------------------------------------------------------------------- 20 IF (A .GT. 2.0) GO TO 30 IF (B .GT. 2.0) GO TO 21 BETALN = GAMLN(A) + GAMLN(B) - GSUMLN(A,B) RETURN 21 W = 0.0 IF (B .LT. 8.0) GO TO 40 BETALN = GAMLN(A) + ALGDIV(A,B) RETURN C C REDUCTION OF A WHEN B .LE. 1000 C 30 IF (B .GT. 1000.0) GO TO 50 N = A - 1.0 W = 1.0 DO 31 I = 1,N A = A - 1.0 H = A/B W = W * (H/(1.0 + H)) 31 CONTINUE W = ALOG(W) IF (B .LT. 8.0) GO TO 40 BETALN = W + GAMLN(A) + ALGDIV(A,B) RETURN C C REDUCTION OF B WHEN B .LT. 8 C 40 N = B - 1.0 Z = 1.0 DO 41 I = 1,N B = B - 1.0 Z = Z * (B/(A + B)) 41 CONTINUE BETALN = W + ALOG(Z) + (GAMLN(A) + (GAMLN(B) - GSUMLN(A,B))) RETURN C C REDUCTION OF A WHEN B .GT. 1000 C 50 N = A - 1.0 W = 1.0 DO 51 I = 1,N A = A - 1.0 W = W * (A/(1.0 + A/B)) 51 CONTINUE BETALN = (ALOG(W) - N*ALOG(B)) + (GAMLN(A) + ALGDIV(A,B)) RETURN C----------------------------------------------------------------------- C PROCEDURE WHEN A .GE. 8 C----------------------------------------------------------------------- 60 W = BCORR(A,B) H = A/B C = H/(1.0 + H) U = -(A - 0.5)*ALOG(C) V = B*ALNREL(H) IF (U .LE. V) GO TO 61 BETALN = (((-0.5*ALOG(B) + E) + W) - V) - U RETURN 61 BETALN = (((-0.5*ALOG(B) + E) + W) - U) - V RETURN END REAL FUNCTION GSUMLN (A, B) C----------------------------------------------------------------------- C EVALUATION OF THE FUNCTION LN(GAMMA(A + B)) C FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2 C----------------------------------------------------------------------- X = DBLE(A) + DBLE(B) - 2.D0 IF (X .GT. 0.25) GO TO 10 GSUMLN = GAMLN1(1.0 + X) RETURN 10 IF (X .GT. 1.25) GO TO 20 GSUMLN = GAMLN1(X) + ALNREL(X) RETURN 20 GSUMLN = GAMLN1(X - 1.0) + ALOG(X*(1.0 + X)) RETURN END REAL FUNCTION BCORR (A0, B0) C----------------------------------------------------------------------- C C EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE C LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A). C IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8. C C----------------------------------------------------------------------- DATA C0/.833333333333333E-01/, C1/-.277777777760991E-02/, * C2/.793650666825390E-03/, C3/-.595202931351870E-03/, * C4/.837308034031215E-03/, C5/-.165322962780713E-02/ C------------------------ A = AMIN1(A0, B0) B = AMAX1(A0, B0) C H = A/B C = H/(1.0 + H) X = 1.0/(1.0 + H) X2 = X*X C C SET SN = (1 - X**N)/(1 - X) C S3 = 1.0 + (X + X2) S5 = 1.0 + (X + X2*S3) S7 = 1.0 + (X + X2*S5) S9 = 1.0 + (X + X2*S7) S11 = 1.0 + (X + X2*S9) C C SET W = DEL(B) - DEL(A + B) C T = (1.0/B)**2 W = ((((C5*S11*T + C4*S9)*T + C3*S7)*T + C2*S5)*T + C1*S3)*T + C0 W = W*(C/B) C C COMPUTE DEL(A) + W C T = (1.0/A)**2 BCORR = (((((C5*T + C4)*T + C3)*T + C2)*T + C1)*T + C0)/A + W RETURN END REAL FUNCTION ALGDIV (A, B) C----------------------------------------------------------------------- C C COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8 C C -------- C C IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY C LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X). C C----------------------------------------------------------------------- DATA C0/.833333333333333E-01/, C1/-.277777777760991E-02/, * C2/.793650666825390E-03/, C3/-.595202931351870E-03/, * C4/.837308034031215E-03/, C5/-.165322962780713E-02/ C------------------------ IF (A .LE. B) GO TO 10 H = B/A C = 1.0/(1.0 + H) X = H/(1.0 + H) D = A + (B - 0.5) GO TO 20 10 H = A/B C = H/(1.0 + H) X = 1.0/(1.0 + H) D = B + (A - 0.5) C C SET SN = (1 - X**N)/(1 - X) C 20 X2 = X*X S3 = 1.0 + (X + X2) S5 = 1.0 + (X + X2*S3) S7 = 1.0 + (X + X2*S5) S9 = 1.0 + (X + X2*S7) S11 = 1.0 + (X + X2*S9) C C SET W = DEL(B) - DEL(A + B) C T = (1.0/B)**2 W = ((((C5*S11*T + C4*S9)*T + C3*S7)*T + C2*S5)*T + C1*S3)*T + C0 W = W*(C/B) C C COMBINE THE RESULTS C U = D*ALNREL(A/B) V = A*(ALOG(B) - 1.0) IF (U .LE. V) GO TO 30 ALGDIV = (W - V) - U RETURN 30 ALGDIV = (W - U) - V RETURN END