c general purpose Bayesian calculation program c c input: poster = fortran function (at the end of this program) to c compute posterior distribution at a given set of parameter c values. First call is to read data and create crude c parameter estimates c data file : first line contains: nseq=# of parallel metrop sequences c ndraws=length of sequences c seed = random seed (0 to 2147483647) c subsequent lines contain data to be read by poster c c output: crude parameter estimates c posterior mode (and approx s.e. from 2nd derivative matrix) c metropolis acceptance rate c metropolis/Gelman-Rubin style output c c parameter statement sets maximum allowed: c number of parameters c number of Metrop sequences c number of draws per sequence (must be even #) c number of values for final inferences = mxseq*mxdraw/2 implicit integer (a-z) parameter(mxparm=50,mxseq=20,mxdraw=100000,mxtot=1000000) double precision parm(mxparm),secder(mxparm,mxparm) double precision invsecdr(mxparm,mxparm),covar(mxparm,mxparm) double precision poster,lh,draws(mxparm+1,mxdraw,mxseq) double precision summ(mxparm+1,10),crudsumm(mxparm+1,3) double precision start(mxseq,mxparm),xp1,xp2,xv c read main control parameters and write first line of output read (5,*) nseq,ndraws,seed,alt,xp1,xp2,xv c read nseq, ndraws, seed, alt, xv=variance for fixed parameters read (4,*) nseq,ndraws,seed,alt,xv if (nseq.gt.mxseq.or.ndraws.gt.mxdraw) then write (6,189) 189 format('too many sequences (>10) or draws (>2000) specified') go to 201 end if nd = ndraws - (ndraws/2) write (6,191) nseq,ndraws,seed,nd 191 format('METROPOLIS PROGRAM',5x,i2,' sequences of length ',i6, . ' with random seed ',i10,/,'metropolis results use last ', . i6,' values in each sequence') write (6,187) xp1,xp2,xv 187 format('mode log(var2): ',f4.1,2x,' mode z(r): ',f4.1,2x . ' var of these: ',f4.2) write (6,192) 192 format(i10) c generate crude estimates and evaluate likelihood at crude estimates nparm = -1 lh = poster(nparm,parm) if (nparm.gt.mxparm) then write (6,188) 188 format('posterior subroutine used too many parameters (>200)') go to 201 end if c save crude estimates for output do 10 i = 1,nparm crudsumm(i,1) = parm(i) 10 continue crudsumm(nparm+1,1) = lh c Newton-Raphson w/ steepest ascent to reach mode c or alternatively enter mode and covariance matrix for t-approx if (alt.eq.0) call getmode(nparm,parm,secder,lh,err) if (alt.eq.1) then read (5,*) (parm(i),i=1,47) parm(48) = xp1 parm(49) = xp2 lh = poster(nparm,parm) read (5,*) ((secder(i,j),j=1,47),i=1,47) do 1605 i = 1,47 do 1605 j = 48,49 secder(i,j) = 0.d0 secder(j,i) = 0.d0 1605 continue secder(48,48) = 1.d0 / xv secder(48,49) = 0.d0 secder(49,48) = 0.d0 secder(49,49) = 1.d0 / xv go to 109 end if if (err.eq.1) then write (6,101) 101 format ('program stopping -- failed to find mode') go to 201 end if c at mode: find second deriv matrix, then cholesky decomp of secder matrix, c (we find upper triangle R with R^T R = second derivative matrix) c then invert cholesky triangle, and finally compute covar matrix c (we find upper triangle R with R^T R = second derivative matrix) c mode (inverse chol x inverse chol^T) call seconder(nparm,parm,secder) 109 call cholesky(secder,nparm,err) if (err.eq.1) then write (6,102) 102 format ('program stopping -- 2nd deriv matrix is not < 0') go to 201 end if call invuptri(secder,nparm,invsecdr) call getcovar(invsecdr,nparm,covar) c collect post mode and approx s.e. for output do 110 i = 1,nparm crudsumm(i,2) = parm(i) crudsumm(i,3) = dsqrt(covar(i,i)) 110 continue crudsumm(nparm+1,2) = lh c c METROPOLIS SECTION OF PROGRAM c includes routines for getting starting values and convergence check c c get starting values for the metropolis sequences call cholesky(covar,nparm,err) call getstart(start,seed,nparm,parm,secder,covar,nseq) c carry out metropolis -- check that acceptance rate is in correct range call metrop(nparm,parm,covar,nseq,start,ndraws,draws,err,seed) if (err.eq.2) write (6,104) 104 format ('never got accept rate in good range in 4 tries') c carry out convergence check on the parallel metropolis sequences call convchk(draws,nseq,ndraws,nparm,summ,err) if (err.eq.3) write (6,103) 103 format ('metrop did not converge (R>1.1 for a param) -- rerun', . ' with more draws') c c write output write (6,197) 197 format(9x,' crud estm',' post mode',' asympt se',10x, . ' metrop t-approxim: mean+95%ci') do 190 i = 1,nparm write (6,198) i,(crudsumm(i,j),j=1,3),(summ(i,j),j=3,5) 198 format ('parm',i3,2x,3(1x,f9.5),10x,3(1x,f9.5)) 190 continue write (6,195) (crudsumm(nparm+1,j),j=1,2),(summ(nparm+1,j),j=3,5) 195 format ('likhood',2x,2(1x,e11.5),20x,3(1x,e11.5)) write (6,192) write (6,196) 196 format(9x,' estim. R ','97.5%ile-R',' 2.5%ile ',' 25%ile ', . ' 50%ile ',' 75%ile ',' 97.5%ile ') do 200 i = 1,nparm write (6,199) i,(summ(i,j),j=1,2),(summ(i,j),j=6,10) 199 format ('parm',i3,2x,7(1x,f9.5)) 200 continue write (6,194) (summ(nparm+1,j),j=1,2),(summ(nparm+1,j),j=6,10) 194 format ('likhood',2x,7(1x,e11.5)) write (6,192) c do 300 i = 1,nseq c do 300 j = 1,ndraws cwrite (6,292) (draws(k,j,i),k=1,nparm+1) c92 format(3f14.8,1x,e14.8) c00 continue 201 stop end subroutine metrop(nparm,parm,covar,nseq,start,ndraws,draws, . err,seed) c generate nseq metropolis chains each of length ndraws implicit integer(a-z) parameter(mxparm=50,mxseq=20,mxdraw=100000) double precision parm(mxparm),covar(mxparm,mxparm) double precision poster,draws(mxparm+1,mxdraw,mxseq) double precision start(mxseq,mxparm),restart(mxseq,mxparm) double precision curr(mxparm),cand(mxparm) double precision curpst,newpst,ratio,chkpost(mxseq) double precision count(mxseq),scale,totchk,unif,rangen err = 0 c Gelman,Gilks,James result suggests this scale for normal scale = 2.4d0/dsqrt(dble(nparm)) c monitor metropolis acceptance rate after ncheck draws (unless ndraws is small) ncheck = 100 if (ndraws.le.2*ncheck) ncheck = ndraws c c metropolis loop -- 1 until time to check on accept rate fail = 0 11 totchk = 0.d0 do 200 i = 1,nseq count(i) = 0 do 1 ii = 1,nparm curr(ii) = start(i,ii) 1 continue curpst = poster(nparm,curr) do 100 j = 1,ncheck c generate candidate from mvn with chosen scale call getcand(nparm,curr,covar,cand,scale,seed) newpst = poster(nparm,cand) ratio = dexp(newpst-curpst) c accept/reject candidate unif = 0.d0 if (ratio.lt.1.d0) unif=rangen(seed) if (unif.lt.ratio) then curpst = newpst count(i) = count(i) + 1 do 10 jj = 1,nparm curr(jj) = cand(jj) 10 continue end if c save draw do 20 jj = 1,nparm draws(jj,j,i) = curr(jj) 20 continue draws(nparm+1,j,i) = curpst 100 continue c if accept rate check is being performed save last draw and posterior if (ncheck.lt.ndraws) then chkpost(i) = curpst do 30 jj = 1,nparm restart(i,jj) = curr(jj) 30 continue totchk = totchk + count(i) end if 200 continue c compute, evaluate, and write metropolis acceptance rate if (ncheck.lt.ndraws) then totchk = totchk / (nseq*ncheck) write (6,250) totchk,ncheck 250 format('metrop accept rate is ',f6.3,' after ',i4,' draws/seq') if (totchk.lt.0.1d0.or.totchk.gt.0.5d0) then fail = fail + 1 if (fail.eq.4) then err = 2 end if if (totchk.gt.0.5d0) scale = scale * 1.25d0 if (totchk.lt.0.1d0) scale = scale / 1.25d0 go to 11 end if c resume metropolis after check (whether accept rate is OK or not) 251 do 290 i = 1,nseq curpst = chkpost(i) do 261 ii = 1,nparm curr(ii) = restart(i,ii) 261 continue do 280 j = ncheck+1,ndraws call getcand(nparm,curr,covar,cand,scale,seed) newpst = poster(nparm,cand) ratio = dexp(newpst-curpst) unif = 0.d0 if (ratio.lt.1.d0) unif=rangen(seed) if (unif.lt.ratio) then curpst = newpst count(i) = count(i) + 1 do 275 jj = 1,nparm curr(jj) = cand(jj) 275 continue end if do 279 jj = 1,nparm draws(jj,j,i) = curr(jj) 279 continue draws(nparm+1,j,i) = curpst 280 continue write (6,285) i 285 format ('done with sequence ',i2) 290 continue end if 291 totchk = 0.d0 do 292 i = 1,nseq totchk = totchk + count(i) 292 continue totchk = totchk / (nseq*ndraws) write (6,250) totchk,ndraws return end subroutine convchk(draws,nseq,ndraws,nparm,summ,err) c assess convergence of metropolis draws implicit integer(a-z) parameter(mxparm=50,mxseq=20,mxdraw=100000,mxtot=1000000) double precision draws(mxparm+1,mxdraw,mxseq) double precision summ(mxparm+1,10),sum,sum2,x,mn,s2,vec(mxtot) double precision summn,sums2,summnsq,sums2sq,desc(5) double precision summns2,summn2s2,temp,fq,qf,tq,qt,df1 double precision W,muhat,B,varW,varB,covWB,postvar double precision sig2hat,varpostvar,postdf,varlodf ni = ndraws/2 nd = ndraws - ni c for each scaler parameter do 1000 ip = 1,nparm+1 count = 0 summn = 0.d0 sums2 = 0.d0 summnsq = 0.d0 sums2sq = 0.d0 summns2 = 0.d0 summn2s2 = 0.d0 c accumulate mean and s.d within each sequence do 100 i = 1,nseq sum = 0.d0 sum2 = 0.d0 do 90 j = 1,nd count = count + 1 x = draws(ip,ni+j,i) vec(count) = x sum = sum + x sum2 = sum2 + x*x 90 continue mn = sum /nd s2 = (sum2 - sum*mn) / (nd-1) c accumulate sums, cross products needed for Gelman-Rubin scale reduction factor summn = summn + mn summnsq = summnsq + mn*mn sums2 = sums2 + s2 sums2sq = sums2sq + s2*s2 summns2 = summns2 + mn*s2 summn2s2 = summn2s2 + mn*mn*s2 100 continue c G-R calculations W = sums2 / nseq muhat = summn / nseq B = nd*(summnsq - summn*muhat)/(nseq-1) varW = (sums2sq - sums2*W)/((nseq-1)*nseq) varB = 2.d0*B*B/(nseq-1) covWB = nd*( ((summn2s2 - summnsq*W)/(nseq-1)) - . ((summns2 - summn*W)*2.d0*muhat/(nseq-1)) )/nseq sig2hat = ((nd-1.d0)*W + B)/nd postvar = sig2hat + B/(nd*nseq) varpostvar = (((nd-1)**2)*varW + ((nseq + 1)**2)*varB/(nseq*nseq) . + 2*(nd-1)*(nseq+1)*covWB/nseq)/(nd*nd) postdf = 2*postvar*postvar / varpostvar tq = qt(postdf) c compute conservative t approximation summ(ip,3) = muhat - tq*sqrt(postvar) summ(ip,4) = muhat summ(ip,5) = muhat + tq*sqrt(postvar) c compute quantiles call analy(vec,count,desc) do 987 jj = 1,5 summ(ip,jj+5) = desc(jj) 987 continue c compute scale reduction factor and .975 point c currently uses a fairly lame F-distn quantile approx varlodf = 2*W*W /varW df1 = nseq - 1.d0 fq = qf(df1,varlodf) summ(ip,1) = dsqrt(postvar*postdf/(W*(postdf-2.d0))) temp = (nseq + 1.d0)*(B/W)/(nseq*nd) temp = ((nd-1.d0)/nd) + temp*fq summ(ip,2) = dsqrt( temp*postdf/(postdf-2.d0) ) if (summ(ip,1).gt.1.1d0) err = 3 1000 continue return end subroutine getmode(nparm,parm,secder,lh,merr) c find mode of posterior distribution parameter(mxparm=50) implicit integer(a-z) double precision parm(mxparm),secder(mxparm,mxparm) double precision tmpparm(mxparm),parminc(mxparm) double precision der1(mxparm),poster double precision toler,ztol,lh,oldlh double precision sum,t data toler/1.d-12/,ztol/1.d-20/ data mxloop/100/ merr = 0 end = 0 lh = -999999.d99 do 10 loop=1,mxloop oldlh = lh c compute vector of first derivatives call firstder(nparm,parm,der1) c compute negative of second derivative matrix numerically call seconder(nparm,parm,secder) c Newton-Raphson step c --get cholesky decomposition of secder (writing over secder) call cholesky(secder,nparm,err) if (err.eq.0) then c --solve for Newton step unless Cholesky decomp fails do 50 i=1,nparm sum = 0.d0 do 49 kk = 1,i-1 sum = sum + secder(kk,i)*parminc(kk) 49 continue parminc(i) = (der1(i) - sum) / secder(i,i) 50 continue do 60 i=1,nparm ii = nparm + 1 - i sum = 0.d0 do 59 kk = ii+1,nparm sum = sum + secder(ii,kk)*parminc(kk) 59 continue parminc(ii) = (parminc(ii) - sum) / secder(ii,ii) 60 continue t = 1.d0 else c can not do Newton-Raphson step -- use steepest ascent step do 70 i=1,nparm parminc(i) = der1(i) 70 continue t = 0.1d0 end if c check that step increases posterior density or reduce step (mainly for c steep ascent) 2003 do 2004 i = 1,nparm 2004 tmpparm(i) = parm(i) + parminc(i)*t lh = poster(nparm,tmpparm) if (lh.le.oldlh) then t = t / 2.d0 if (t.lt.ztol) then t = 0.d0 go to 2009 end if go to 2003 end if c found a good step 2009 do 2010 i = 1,nparm 2010 parm(i) = parm(i) + t*parminc(i) c test for convergence if (lh.lt.(oldlh+toler)) then end = 1 go to 101 end if 10 continue merr = 1 write (6,102) mxloop,oldlh,lh 102 format(' did not find max in ',i3,' steps',/, . 'poster dens for last two steps:',2e15.8) 101 return end subroutine getstart(start,seed,nparm,parm,secder,covar,nseq) c get starting values for metropolis by impt wtd sampling from t_4 approx parameter(mxparm=50,mxseq=20,mximpt=2000) implicit integer(a-z) double precision parm(mxparm),secder(mxparm,mxparm) double precision covar(mxparm,mxparm) double precision gaudev,gamdev,sum,bigset(mximpt,mxparm) double precision norms(mxparm),poster,new(mxparm),tdens double precision rangen,uni,start(mxseq,mxparm) double precision fac,sumwt,thewt(mximpt),bigwt(mximpt) double precision aa,bb data four/4/,two/2/ c generate 1000 draws from t_4 approx at post mode sumwt = 0.d0 do 1000 i = 1,mximpt do 10 ii = 1,nparm norms(ii) = gaudev(seed) 10 continue fac = dsqrt(2.d0/gamdev(two,seed)) do 100 ii = 1,nparm sum = 0.d0 do 99 j = 1,ii sum = sum + covar(j,ii)*norms(j) 99 continue new(ii) = parm(ii) + sum*fac bigset(i,ii) = new(ii) 100 continue c compute importance ratio thewt(i) = dexp(poster(nparm,new)+5900- . tdens(new,parm,secder,nparm,four)) sumwt = sumwt + thewt(i) bigwt(i) = sumwt 1000 continue c sample nseq starting values using importance ratios as weights (w/out replac) do 2000 i = 1,nseq uni = rangen(seed)*sumwt c binary search to locate sampled value call locate(bigwt,mximpt,uni,ix) do 1100 j = 1,nparm start(i,j) = bigset(ix,j) 1100 continue sumwt = sumwt - thewt(ix) do 1200 k = ix,mximpt bigwt(k) = bigwt(k) - thewt(ix) 1200 continue 2000 continue return end subroutine getcand(nparm,curr,vcov,cand,scale,seed) c generate candidate value for metropolis -- using normal jumping distn with c s.d. multiplied by scale parameter(mxparm=50) implicit integer(a-z) double precision curr(mxparm),cand(mxparm),vcov(mxparm,mxparm) double precision scale,gaudev,norms(mxparm),sum do 10 i = 1,nparm norms(i) = gaudev(seed) 10 continue do 100 i = 1,nparm sum = 0.d0 do 99 j = 1,i sum = sum + vcov(j,i)*norms(j) 99 continue cand(i) = curr(i) + scale*sum 100 continue return end double precision function tdens(x,mu,sigma,dim,df) c evaluate the t-density used to generate Metropolis starts implicit integer (a-z) parameter(mxparm=50) double precision x(mxparm),mu(mxparm),sigma(mxparm,mxparm) double precision sum,sumvsq,rdf rdf = df sumvsq = 0.d0 do 60 i = 1,dim sum = 0.d0 do 59 kk = i,dim sum = sum + sigma(i,kk)*(x(kk)-mu(kk)) 59 continue sumvsq = sumvsq + sum*sum 60 continue tdens = -0.5d0*(dim+df)*dlog(1.d0 + (sumvsq/rdf)) return end subroutine firstder(nparm,parm,der1) c numerically evaluate vector of first derivatives of posterior distn implicit integer (a-z) parameter(mxparm=50) double precision step,delta,parm(mxparm),der1(mxparm) double precision l1,l2,poster data step/1.d-04/ do 10 i=1,nparm delta = step*parm(i) parm(i) = parm(i) + delta l1 = poster(nparm,parm) parm(i) = parm(i) - 2.d0*delta l2 = poster(nparm,parm) parm(i) = parm(i) + delta der1(i) = (l1 - l2)/(2.d0*delta) 10 continue return end subroutine seconder(nparm,parm,secder) c numerically evaluate matrix of second derivatives of posterior distn c (actually finds negative of second derivs --> inverse of asymp covar matrix) implicit integer (a-z) parameter(mxparm=50) double precision parm(mxparm),secder(mxparm,mxparm) double precision step,delta1,delta2,poster double precision l1,l2,l3,l4 data step/1.d-04/ do 10 i = 1,nparm delta1 = step*parm(i) do 11 j = i,nparm delta2 = step*parm(j) parm(i) = parm(i) + delta1 parm(j) = parm(j) + delta2 l1 = poster(nparm,parm) parm(j) = parm(j) - 2.d0*delta2 l2 = poster(nparm,parm) parm(i) = parm(i) - 2.d0*delta1 l3 = poster(nparm,parm) parm(j) = parm(j) + 2.d0*delta2 l4 = poster(nparm,parm) secder(i,j) = -1.d0*(l1 + l3 - l2 - l4)/(4.d0*delta1*delta2) parm(i) = parm(i) + delta1 parm(j) = parm(j) - delta2 11 continue secder(j,i) = secder(i,j) 10 continue return end subroutine getcovar(invchlpr,n,covar) c compute asymptotic covar matrix of normal/t approx c --- if R is cholesky decomp of second deriv matrix and T = R^(-1) c --- then this finds T T^t implicit integer (a-z) parameter(mxparm=50) double precision invchlpr(mxparm,mxparm),covar(mxparm,mxparm) double precision sum do 200 i = 1,n do 199 j = i,n sum = 0.d0 do 198 k = j,n sum = sum + invchlpr(i,k)*invchlpr(j,k) 198 continue covar(i,j) = sum if (j.ne.i) covar(j,i) = sum 199 continue 200 continue return end subroutine invuptri(x,n,y) c compute inverse of upper triangular matrix (used here to invert cholesky c triangle of second deriv matrix) implicit integer (a-z) parameter(mxparm=50) double precision x(mxparm,mxparm),y(mxparm,mxparm),sum do 200 i = 1,n j = n + 1 - i y(j,j) = 1.d0 / x(j,j) do 200 ii=1,j-1 k = j - ii sum = 0.d0 do 199 iii = k+1,j sum = sum + x(k,iii)*y(iii,j) 199 continue y(k,j) = -1.d0*sum/x(k,k) 200 continue return end subroutine cholesky(x,n,err) c compute cholesky decomposition of the matrix x implicit integer(a-z) parameter(mxparm=50) double precision x(mxparm,mxparm),sum err = 0 do 32 i = 1,n sum = 0.d0 do 29 kk = 1,i-1 sum = sum + x(kk,i)*x(kk,i) 29 continue if (x(i,i).le.sum) then err = 1 return end if x(i,i) =dsqrt(x(i,i)-sum) do 31 j = i+1,n sum = 0.d0 do 30 kk = 1,i-1 sum = sum + x(kk,j)*x(kk,i) 30 continue x(i,j) = (x(i,j) - sum)/x(i,i) 31 continue 32 continue do 40 i = 2,n do 40 j = 1,i-1 x(i,j) = 0.d0 40 continue return end subroutine locate(xx,n,x,j) c binary search routine (used to carry out importance weighted samples) parameter(mximpt=2000) implicit integer (a-z) double precision xx(mximpt),x jl = 0 ju = n 10 if ((ju-jl).gt.1.) then jm = (ju + jl)/2 if (x.gt.xx(jm)) then jl = jm else ju = jm end if go to 10 end if j = ju return end subroutine analy(x,n,desc) c sorts vector and obtains quantiles (2.5, 25, 50, 75, 97.5) parameter(mxtot=1000000) implicit integer(a-z) double precision x(mxtot) double precision desc(5),ixr,r call hpsort(n,x) ix = (n+1)/40 ixr = dble(n+1)/40 r = ixr - ix desc(1) = (1.0-r)*x(ix) + r*x(ix+1) ix = (n+1)/4 ixr = dble(n+1)/4 r = ixr - ix desc(2) = (1.0-r)*x(ix) + r*x(ix+1) ix = (n+1)/2 ixr = dble(n+1)/2 r = ixr - ix desc(3) = (1.0-r)*x(ix) + r*x(ix+1) ix = 3*(n+1)/4 ixr = 3*dble(n+1)/4 r = ixr - ix desc(4) = (1.0-r)*x(ix) + r*x(ix+1) ix = 39*(n+1)/40 ixr = 39*dble(n+1)/40 r = ixr - ix desc(5) = (1.0-r)*x(ix) + r*x(ix+1) return end subroutine hpsort(n,ra) c heap sort from Numerical recipes implicit integer (a-z) parameter(mxtot=1000000) double precision rra,ra(mxtot) l = (n/2) + 1 ir = n 10 continue if (l.gt.1) then l = l - 1 rra=ra(l) else rra=ra(ir) ra(ir)=ra(1) ir = ir - 1 if (ir.eq.1) then ra(1) = rra return end if end if i = l j = l + l 20 if (j.le.ir) then if (j.lt.ir) then if (ra(j).lt.ra(j+1)) j = j + 1 end if if (rra.lt.ra(j)) then ra(i)=ra(j) i = j j = j + j else j = ir + 1 end if go to 20 end if ra(i) = rra go to 10 end double precision function rangen(ix) c portable uniform random number generator -- mult=16807, base=2**31 double precision m integer a,p,ix,b15,b16,xhi,xalo,leftlo,fhi,k data a/16807/,b15/32768/,b16/65536/,p/2147483647/ data m/4.6566128752458e-10/ xhi = ix/b16 xalo = (ix-xhi*b16)*a leftlo = xalo/b16 fhi = xhi*a + leftlo k = fhi / b15 oldix = ix ix = (((xalo - leftlo*b16) - p) + (fhi-k*b15)*b16) + k if (ix.lt.0) ix = ix + p rangen = dble(ix)*m if (rangen.ge.1.d0.or.rangen.le.0.d0) write (6,1) 1 format('apparent problem in rangen') return end double precision function gamdev(alph,iseed) c generate gamma random variables (shape parameter = alph is an integer) c from Numerical Recipes implicit integer (a-z) double precision x,rangen,y,s,e,pi,alphm1 data pi/3.141592654d0/ if (alph.le.0) then write (6,1) 1 format ('problem in gamdev') x = 0.d0 else if (alph.lt.6) then x = 1.d0 do 11 j = 1,alph x = x*rangen(iseed) 11 continue x = -dlog(x) else 21 y = tan(pi*rangen(iseed)) alphm1 = alph - 1.d0 s = dsqrt(2.d0*alphm1+1.d0) x = s*y + alphm1 if (x.le.0.d0) go to 21 e = (1.d0+y*y)*dexp(alphm1*dlog(x/alphm1)-s*y) if (rangen(iseed).gt.e) go to 21 end if gamdev = x return end double precision function gaudev(ix) c generate normal random variable implicit integer (a-z) double precision v1,v2,r,gset,fac,rangen data iset/0/ if (iset.eq.0) then 1 v1 = 2.d0*rangen(ix) - 1.d0 v2 = 2.d0*rangen(ix) - 1.d0 r = v1**2 + v2**2 if (r.ge.1.d0) go to 1 fac = dsqrt(-2.d0*dlog(r)/r) gset = v1*fac gaudev= v2*fac iset = 1 else gaudev = gset iset = 0 endif return end double precision function qt(df) c from Abramowitz and Stegun -- inverse function for t-distn c currently hardwired for .975 pct point -- to change just replace xp implicit integer (a-z) double precision df,xp,xp3,xp5,xp7,xp9 data xp/1.959964d0/ xp3 = xp*xp*xp xp5 = xp3*xp*xp xp7 = xp5*xp*xp xp9 = xp7*xp*xp qt = xp + (xp3+xp)/(4.d0*df) + . (5*xp5 + 16*xp3 + 3*xp)/(96.d0*df*df) + . (3*xp7 + 19*xp5 + 17*xp3 - 15*xp)/(384.d0*df*df*df) + . (79*xp9 + 776*xp7 + 1482*xp5 - 1920*xp3 - 945*xp)/ . (92160.d0*df*df*df*df) return end double precision function qf(df1,df2) c from Abramowitz and Stegun -- approx to inverse function for F-distn c approximation doesn't look too good c currently hardwired for .975 pct point -- to change just replace xp double precision df1,df2,a,b,lam,h,w,xp data xp/1.959964d0/ a = df2 / 2.d0 b = df1 / 2.d0 lam = (xp*xp - 3.d0) / 6.d0 h = 2.d0 / ( (1.d0/(df1-1.d0)) + (1.d0/(df2-1.d0)) ) w = xp*sqrt(h + lam)/h - (lam+(5.d0/6.d0)-(2.d0/(3.d0*h)))* . ((1.d0/(df1-1.d0))-(1.d0/(df2-1.d0))) qf = dexp(2.d0*w) return end double precision function invlogt(x) implicit integer (a-z) double precision x,t t= dexp(x) invlogt = t/(1.d0+t) return end double precision function logit(y,n) implicit integer (a-z) double precision x,yy yy = y if (y.eq.0) yy = 0.5d0 if (y.eq.n) yy = n - 0.5d0 x = yy/n logit = dlog(x/(1.d0-x)) return end