subroutine antith(kmax,k,irsimBIG,it,v0,a) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c ANTITH is a routine that creates a "grid" of antithetic points S on the c unit sphere, starting from a random basis provided by RNDBASE. c These routines will ordinarily be called for each observation. To c avoid "chatter" when parameter values are changed, the seed for the c random number generator for each observation should be stored so that c the same outputs can be regenerated in iteration to parameter estimates. c ccc k=m c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension v0(kmax,irsimBIG),a(kmax,2*kmax+4) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c do 10 i = 1,k do 10 j = 1,k v0ij = v0(i,j) a(i,j) = v0ij a(i,k+j) = -v0ij 10 continue if (it .le. 1) then go to 999 else xit=it xk=k j = 2 + intfloor((0.5d0*xit-xk)/(2.d0*xk*(xk-1.d0))) dj = j do 100 i = 1,k-1 do 200 ii = i+1,k jj = 1 300 continue djj = jj dj = djj / dj if (jj.lt.j) then do 400 n = 1,k v1n = v0(n,i)*dcos(1.5707963d0*dj) v2n = v0(n,ii)*dsin(1.5707963d0*(1.0d0-dj)) a(n,k+1) = v1n + v2n a(n,k+2) = v1n - v2n a(n,k+3) = -v1n + v2n a(n,k+4) = -v1n - v2n 400 continue jj = jj + 1 go to 300 endif 200 continue 100 continue endif 999 continue return end subroutine arse(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c P. Acceptance/Rejection Simulator (ARS) c The ARS simulator is given in a version using the truncated c bilateral exponential comparison distribution (ARSE), and in a c version using a recursive truncated comparison distribution (ARSR). c PARM[.,1] is assumed to contain seeds for the uniform random number c generator RNDUS called by the routines; there is one seed for each c repetition R. PARM[R+1,1] is assumed to contain a censoring limit c for rejection. Should only compute HC. c c The ARS method provides a mechanism for drawing from a conditional c density when transformations from uniform or standard normal variates c are not available. We want to sample from f(x/a). We generate x c from g(x). In our case g(x) = ARSE. Generate tsi from uniform. c Generate unconditional f(x) and bound alpha. We accept x if c tsi < f(x)/(g(x)*alpha) c x has conditional density f(x/a) = HC. We set x = VD. c c NB: Upon exit, the first element of the BLANK common block, contains c -- the proportion of non-censored draws c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / e(mmaxloc,irmaxloc),vd(mmaxloc), / bxmu(mmaxloc),axmu(mmaxloc),temp(mmaxloc),d(mmaxloc), / pa(mmaxloc),pb(mmaxloc),xmatmm(mmaxloc,mmaxloc), / tempt(1,mmaxloc) common/times/hseclim c mp1=m+1 r=ir mmaxp1 = mmax + 1 rscale = 1.0d0/r q = 0.0d0 do 100 j = 1,m do 101 i = 1,m xmatmm(i,j) = wi(i,j) 101 continue 100 continue call invv(xmatmm,xmatmm,mmaxloc,mmaxloc,vd,m,npr,nzr,nnr,rcond, * *999) c ccc-----note:vd contains the eigen values of wi in descending order----- c call sorts(vd,m) alp = 0.0d0 do 200 i = 1,m wii = w(i,i) si = dsqrt(wii+1.d-100) xmui = xmu(i) bxmui = b(i) - xmui axmui = a(i) - xmui tempi = -3.0d0 * si di = 3.0d0 * si t1 = 0.0d0 t2 = 0.0d0 t3 = 0.0d0 t4 = 0.0d0 if (bxmui.lt.tempi) t1=1.0d0 if (axmui.gt.di) t2=1.0d0 if (bxmui.lt.di) t3=1.0d0 if (axmui.gt.tempi) t4=1.0d0 t1 = t1*(bxmui-2.0d0*si) t1 = t1*(bxmui-2.0d0*si) + t2*axmui * + (1.0d0-t1-t2)*(tempi+t4*(axmui+di)) t2 = t1*bxmui + t2*(axmui+2.0d0*si) * + (1.0d0-t1-t2)*(di+t3*(bxmui+tempi)) bb = 0.0d0 aa = 0.0d0 if (t2.lt.0.0d0) aa=1.0d0 if (t1.gt.0.0d0) bb=1.0d0 si = 1.0d0 - aa - bb tempi = (t2-t1)*(0.5d0-si/4.0d0)/(wii+1.d-100) alp = alp + tempi**2.0d0 bxmu(i) = bxmui axmu(i) = axmui temp(i) = tempi 200 continue alp = alp/(2.0d0*vd(1)+1.d-100) call bexp(mmaxloc,m,axmu,temp,pa) call bexp(mmaxloc,m,bxmu,temp,pb) jl = parm(ir+1,1) do 400 j = 1,mp1 do 401 i = 1,mp1 hc(i,j) = 0.0d0 401 continue 400 continue timelim = hsec1990() do 500 i = 1,ir isd = idint(parm(i,1)) do 600 j = 1,jl do 610 k = 1,m cvahfix uk = unifrm(isd) uk = randu(isd) d(k) = uk*pa(k) + (1.0d0-uk)*pb(k) 610 continue call bexpi(mmaxloc,m,1,d,temp,pa,pb,vd,d) cvahfix uzscal = unifrm(isd) uzscal = randu(isd) call matptv(vd,mmaxloc,wi,mmax,tempt,1,m,1,m) scalz1 = 0.0d0 scalz2 = 0.0d0 do 620 ii = 1,m vdii = vd(ii) scalz1 = scalz1 + tempt(1,ii)*vdii scalz2 = scalz2 + temp(ii)*dabs(vdii) 620 continue scale = -(scalz1/2.0d0) + scalz2 - alp if ( dlog(uzscal) .lt. scale ) then q = q + rscale jl = j go to 88 endif 600 continue 88 continue call matpv(wi,mmax,vd,mmaxloc,vd,mmaxloc,m,m,1) do 700 j = 1,m tempt(1,j) = vd(j) 700 continue call matpv(vd,mmaxloc,tempt,1,xmatmm,mmaxloc,m,1,m) hc(1,1) = hc(1,1)+1.d0 do 800 j=1,m vdj = vd(j) hc(1,j+1) = hc(1,j+1)+vdj hc(j+1,1) = hc(j+1,1)+vdj do 801k=1,m hc(k+1,j+1) = hc(k+1,j+1)+(xmatmm(k,j) - wi(k,j))*0.5d0 801 continue 800 continue if ( (hsec1990()-timelim) .gt. hseclim ) then c give up r = i rscale = 1.0d0/r go to 99 endif 500 continue 99 continue do 910 j = 1,m+1 do 911 i = 1,m+1 hc(i,j) = hc(i,j)*rscale hu(i,j) = -999.0d0 911 continue 910 continue p(1) = -999.0d0 ccc print *, 'arse proportion non-censored', q e(1,1)=q go to 888 999 continue print *, 'error in routine finding eigen values' 888 continue return end subroutine arsr(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c see ARSE for documentation. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / uu(mmaxloc+1,irmaxloc),tt(mmaxloc), / wic(mmaxloc,mmaxloc),v(mmaxloc),vt(1,mmaxloc), / htmm(mmaxloc,mmaxloc) common/times/hseclim c mp1 = m + 1 mmaxp1 = mmax + 1 mmaxlcp1 = mmaxloc + 1 r=ir rscale = 1.0d0/r do 10 j = 1,m do 11 i = 1,m hc(i,j) = 0.0d0 11 continue 10 continue jl = parm(ir+1,1) chw11 = chw(1,1) xmu1 = xmu(1) q = 0.0d0 timelim = hsec1990() do 100 i = 1,ir isd = idint(parm(i,1)) do 200 j = 1,jl ta = pheev( (a(1)-xmu1) / (chw11+1.d-100) ) tb = pheev( (b(1)-xmu1) / (chw11+1.d-100) ) call mrunif(uu,mmaxlcp1,mp1,jl,isd) do 210 jk = 1,jl do 211 ik = 1,mp1 cvahfix uu(i,j) = unifrm(isd) uu(i,j) = randu(isd) 211 continue 210 continue tt(1) = phinv( uu(1,j)*ta + (1.0d0 - uu(1,j))*tb ) wgt = tb - ta alp = wgt do 300 k = 2,m km1 = k - 1 ak = a(k) bk = b(k) xmuk = xmu(k) chwkk = chw(k,k) uukj = uu(k,j) chw1tt = 0.0d0 do 301 jk = 1,km1 chw1tt = chw1tt + chw(k,jk)*tt(jk) 301 continue ta = pheev( (ak-xmuk-chw1tt) / (chwkk+1.d-100) ) tb = pheev( (bk-xmuk-chw1tt) / (chwkk+1.d-100) ) tt(k) = phinv(uukj*ta+(1.0d0-uukj)*tb) wgt = wgt*(tb-ta) cdfn = pheev((bk-ak)/(2.0d0*chwkk+1.d-100)) alp = alp*(2.0d0*cdfn - 1.0d0) 300 continue if ( uu(mp1,j) .lt. wgt/alp ) then q = q + rscale jl = j go to 310 endif 200 continue 310 continue call matpv(wi,mmax,chw,mmax,wic,mmaxloc,m,m,m) call matpv(wic,mmaxloc,tt,mmaxloc,v,mmaxloc,m,m,1) do 400 ik = 1,m vt(1,ik) = v(ik) 400 continue call matpv(v,mmaxloc,vt,1,htmm,mmaxloc,m,1,m) hc(1,1) = hc(1,1)+1.d0 do 500 j=1,m vj = v(j) hc(1,j+1) = hc(1,j+1)+vj hc(j+1,1) = hc(j+1,1)+vj do 501 k=1,m hc(k+1,j+1) = hc(k+1,j+1)+(htmm(k,j) - wi(k,j))*0.5d0 501 continue 500 continue if ( (hsec1990()-timelim) .gt. hseclim ) then c give up r=i rscale = 1.0d0/r goto 900 endif 100 continue 900 continue do 910 j = 1,m+1 do 911 i = 1,m+1 hc(i,j) = hc(i,j)*rscale hu(i,j) = -999.0d0 911 continue 910 continue p(1) = -999.0d0 ccc print *, 'arsr proportion non-censored', q uu(1,1)=q return end subroutine aus(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c S. Asymptotically Unbiased Simulator (AUS) c The AUS simulator is obtained as the product of an unbiased c simulator of H and an asymptotically unbiased simulator of 1/P. The c proc is defined with the GHK simulator used for H and for the initial c independent simulations of P(B)/P(A), where A is the event that the c first component is in the rectangle bounds. The proc assumes that c PARM(.,1) contains a vector of R seeds for the normal random number c generator RNDNS, PARM(1,2) contains the number of initial GHK c simulators of P used in the simulation of 1/P, PARM(2,2) contains the c maximum number of trials using the crude frequency simulator that are c made before the process is censored. If PARM(2,2) is made large (say c 1.e50), then this is computationally the same as the Sequential c Unbiased Simulator (SUS). PARM(3,2) contains a seed for RNDUS. c Gives HU and HC. c c This does continuous parameter estimation of 1/p. c c NB: Upon exit, the first element of the BLANK common block, contains c -- the maximum number of CFS draws allowed in AUS. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,ir),parm(irmax+1,2), / hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / f(mmaxloc-1),g(mmaxloc-1,mmaxloc-1),vt(mmaxloc-1), / a2(mmaxloc-1),b2(mmaxloc-1) c common/times/hseclim c mm1 = m - 1 mmaxlcm1 = mmaxloc - 1 mp1 = m + 1 chw11 = chw(1,1) xmu1 = xmu(1) do 101 j = 2,m f(j-1) = w(j,1)/(w(1,1)+1.d-100) do 102 i = 2,m g(i-1,j-1) = chw(i,j) 102 continue 101 continue ta = pheev((a(1)-xmu1)/(chw11+1.d-100)) tb = pheev((b(1)-xmu1)/(chw11+1.d-100)) jlim = parm(2,2) ccc print *, 'in Aus: max. num. of cfs draws = ', jlim isd = idint(parm(3,2)) ii = 2 call fghk(ii,mmax,m,xmu,w,wi,chw,a,b,ir,isd,pa,pt,hu) pinv = 1.0d0 q = 1.0d0 ii = 1 do 20 j = 1,parm(1,2) call fghk(ii,mmax,m,xmu,w,wi,chw,a,b,ir,isd,pa,rr,hc) q = q*(1.0d0-pt/(pa+1.d-100)) pinv = pinv + q 20 continue ppinv = 0.0d0 timelim = hsec1990() do 30 i = 1,ir rr = 0.0d0 isdn = idint(parm(i,1)) do 40 j = 1,jlim cvahfix uu = unifrm(isd) uu = randu(isd) tt = chw11*phinv(uu*ta + (1.0d0-uu)*tb) do 49 k=1,mm1 a2(k)=ranorm(0.d0,1.d0,isdn) 49 continue do 50 k = 1,mm1 temp=0.d0 do 51 kk = 1,mm1 temp=temp+g(k,kk)*a2(kk) 51 continue vt(k)=f(k)*tt+temp 50 continue do 600 k = 1,mm1 temp = vt(k) a2k = a(k+1) - xmu(k+1) b2k = b(k+1) - xmu(k+1) a2(k) = temp - a2k b2(k) = temp - b2k 600 continue a2zmin = fmin(a2,mm1) b2zmax = fmax(b2,mm1) ccc if ( a2zmin .gt. 0.0d0 .and. b2zmax .lt. 0.0d0 ) then if ( a2zmin .gt. 0.0d0 .and. b2zmax .lt. 0.0d0 * .or. (hsec1990()-timelim) .gt. hseclim ) then jlim = j go to 700 else rr = rr + 1.0d0 endif 40 continue 700 continue ppinv = ppinv + pinv + q*rr ccc print *,'ppinv',ppinv if ( (ppinv .gt. 1.d+50) .or. * (hsec1990()-timelim) .gt. hseclim ) then c give up either to prevent overflow or because of time: ir = i go to 900 endif 30 continue 900 continue ccc ppinv = ppinv/(ir*pa+1.d-100) temp=ir if(pa.gt.0.d0) then temp=temp*pa ppinv = ppinv/temp do 910 j = 1,mp1 do 910 i = 1,mp1 hc(i,j) = hu(i,j)*ppinv 910 continue else do 911 j = 1,mp1 do 911 i = 1,mp1 hc(i,j) = -999.d0 911 continue endif p = hu(1,1) c place jlim=parm(2,2) in common (must use double precision): f(1)=parm(2,2) ccc print *,'jlim parm(2,2) f(1)',jlim,parm(2,2),f(1) return end subroutine bexp(mmax,m,v,xlam,p) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c bilateral exponential cdf. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension v(mmax),xlam(mmax),p(mmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / e(mmaxloc,irmaxloc) c do 10 i = 1,m if (v(i) .gt. 0.0d0) then p(i) = 1.0d0 - 0.5d0*dexp(xlam(i)*(-dabs(v(i)))) else p(i) = 0.5d0*dexp(xlam(i)*(-dabs(v(i)))) endif 10 continue return end subroutine bexpi(mmax,m,ir,p,xlam,pa,pb,x,d) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c inverse truncated bilateral exponential cdf. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension p(mmax,ir),xlam(mmax),pa(mmax,ir),pb(mmax,ir), / x(mmax,ir),d(mmax,ir) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // e(mmaxloc,irmaxloc) c do 10 j = 1,ir do 20 i = 1,m pij = p(i,j) xlami = xlam(i) if ( pij .gt. 0.5d0 ) then phalf = 1.0d0 else phalf = 0.0d0 endif pp = pij + phalf*(1.0d0-2.0d0*pij) x(i,j) = (1.0d0-2.0d0*phalf)*dlog(2.0d0*pp+1.d-100) * /(xlami+1.d-100) d(i,j) = xlami*pp /( pb(i,j)-pa(i,j)+1.d-100 ) 20 continue 10 continue return end subroutine cf(irmax,m,ir,alp,bet,gam,y0,y1,y2) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c see PCF for documentation. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension alp(irmax),bet(irmax),gam(irmax),y0(irmax), / y1(irmax),y2(irmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / t(irmaxloc),z0(irmaxloc), / z1(irmaxloc),z2(irmaxloc) c do 100 i = 1,ir alpi = alp(i) + 1.d-100 gami = gam(i) beti = bet(i) dalp = dsqrt(alpi) ti = dexp(-0.5d0*alpi*(gami-beti/alpi)**2.0d0) * /alpi z0i = 2.506628274631 * (pheev(beti/dalp)) / dalp y0i = -2.506628274631 * *(pheev((beti-alpi*gami)/dalp)) / dalp + z0i z1i = (dexp(-beti**2/(2.0d0*alpi))) / alpi y1i = z1i - ti + y0i*beti/alpi z1i = z1i + z0i*beti/alpi z2i = z1i*beti/alpi + z0i/alpi y2i = y1i*beti/alpi + y0i/alpi - gami*ti alp(i) = alpi t(i) = ti y0(i) = y0i y1(i) = y1i y2(i) = y2i z0(i) = z0i z1(i) = z1i z2(i) = z2i 100 continue k = m-1 if (k .gt. 0) then j = 2 200 continue if (j .lt. k+2) then dj = j do 202 i = 1,ir alpi = alp(i) beti = bet(i) gji = gam(i)**j y0i = y1(i) y1i = y2(i) y2i = y1i*beti/alpi + y0i*dj/alpi * - gji*t(i) z0i = z1(i) z1i = z2(i) z2i = z1i*beti/alpi + z0i*dj/alpi y0(i) = y0i y1(i) = y1i y2(i) = y2i z0(i) = z0i z1(i) = z1i z2(i) = z2i 202 continue j = j+1 goto 200 endif endif return end subroutine cfs(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c I. Crude Frequency Simulator (CFS) c This simulator assumes that U is an array of standard normal c variates, independent within each column, and either independent or c antithetic across columns. It does not use PARM. It returns both HU c and HC. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,m),wi(mmax,m),chw(mmax,m), / u(mmax,ir),parm(irmax+1,2), / p(1),hu(mmax+1,m+1),hc(mmax+1,m+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / vd(mmaxloc,irmaxloc), / f(irmaxloc),xmatmir(mmaxloc,irmaxloc), / vdsq(mmaxloc,mmaxloc),vdt(irmaxloc,mmaxloc), / vdd(mmaxloc) c mp1 = m+1 r = ir rscale = 1.0d0/r call matpv(chw,mmax,u,mmax,vd,mmaxloc,m,m,ir) do 10 j=1,ir do 11 i=1,m aa = 0.0d0 bb = 0.0d0 axmu = a(i)-xmu(i) bxmu = b(i)-xmu(i) if (vd(i,j).lt.bxmu) bb = 1.0d0 if (vd(i,j).gt.axmu) aa = 1.0d0 xmatmir(i,j) = aa*bb 11 continue 10 continue sum = 0.0d0 do 20 j = 1,ir fj = fmin(xmatmir(1,j),m) f(j) = fj sum = sum + fj do 21 i = 1,m xmatmir(i,j) = vd(i,j)*fj 21 continue 20 continue p(1) = sum*rscale call matpv(wi,mmax,xmatmir,mmaxloc,vd,mmaxloc,m,m,ir) do 30 j = 1,m sum = 0.0d0 do 31 i = 1,ir temp = vd(j,i) vdt(i,j) = temp sum = sum + temp 31 continue vdd(j) = sum*rscale 30 continue call matpv(vd,mmaxloc,vdt,irmaxloc,vdsq,mmaxloc,m,ir,m) hu(1,1) = p(1) do 40 j=1,m hu(1,j+1) = vdd(j) hu(j+1,1) = vdd(j) do 41 i=1,m hu(i+1,j+1) = (vdsq(i,j)*rscale - p(1)*wi(i,j))*0.5d0 41 continue 40 continue if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 50 j=1,mp1 do 51 i=1,mp1 hc(i,j)=hu(i,j)*pscale 51 continue 50 continue else do 60 j=1,mp1 do 61 i=1,mp1 hc(i,j)=-999.d0 61 continue 60 continue endif return end subroutine dcs(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c O. Deak Chi-Square Simulator (DCS) c c The DCS assumes that U is an array whose columns are points in c the unit sphere, antithetic across columns. c The routine assumes that PARM(1,1) = DET(W). c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / aasp(mmaxloc+1,irmaxloc), / bbsp(mmaxloc,irmaxloc),qu(irmaxloc),ql(irmaxloc), / temp(1,mmaxloc),tempv(mmaxloc),htmm(mmaxloc,mmaxloc), / tempm(mmaxloc,irmaxloc) cvah / tempm(mmaxloc,irmaxloc),det(2),kpvt(mmaxloc) c common/times/hseclim c cvah dimension inert(3) c dm = m mp1 = m + 1 do 100 j = 1,m do 101 i = 1,m tempm(i,j) = w(i,j) 101 continue 100 continue cvah det(w) must be passed as parm(1,1) detw=parm(1,1) cvah call dsifa(tempm,mmaxloc,m,kpvt,info) cvah if(info.eq.0) then cvah info=010 cvah call dsidi(tempm,mmaxloc,m,kpvt,det,inert,tempv,info) cvah else cvah print *,'W matrix singular: in DCS routine. Aborting...' cvah stop cvah endif cvah detw = det(1)*10.0d0**det(2) dk = 1.0d0/dsqrt(detw + 1.d-100) gammo2 = dgamma(dm/2.0d0) + 1.d-100 dk1 = dk * dgamma((dm+1.d0)/2.d0) * dsqrt(2.0d0) / gammo2 dk2 = dk * 2.0d0 * dgamma((dm+2.d0)/2.d0) / gammo2 do 200 j = 1,ir do 201 i = 1,m ai = a(i) bi = b(i) xmui = xmu(i) uij = u(i,j) if ( uij .eq. 0.0d0 ) then uu = uij + 1.d-100 else uu = uij endif sn = 0.0d0 sp = 0.0d0 if (uij .lt. 0.0d0) sn = 1.0d0 if (uij .gt. 0.0d0) sp = 1.0d0 aasp(i,j) = ((ai-xmui)*sp+(bi-xmui)*sn) / (uu+1.d-100) bbsp(i,j) = ((bi-xmui)*sp+(ai-xmui)*sn) / (uu+1.d-100) 201 continue aasp(m+1,j) = 0.0d0 qlj = fmax(aasp(1,j),mp1) quj = fmin(bbsp(1,j),m) if (qlj.gt.quj) then qu(j) = qlj else qu(j) = quj endif ql(j) = qlj 200 continue timelim = hsec1990() dfm = m dfmp1 = mp1 dfmp2 = m+2 do 500 j = 1,mp1 do 501 i = 1,mp1 hu(i,j) = 0.0d0 501 continue 500 continue do 600 j = 1,ir call matptv(u(1,j),mmax,wi,mmax,temp,1,m,1,m) t = 0.0d0 do 610 i = 1,m t = t + temp(1,i)*u(i,j) 610 continue tr = 1.0d0/(dsqrt(t)+1.d-100) tt = tr**m if ( qu(j) .gt. ql(j) ) then xzql = t*ql(j)**2 xzqu = t*qu(j)**2 d0zscale = - chicdf(xzql,dfm,ier) + chicdf(xzqu,dfm,ier) d0 = dk*tt*d0zscale d1zscale = - chicdf(xzql,dfmp1,ier) + chicdf(xzqu,dfmp1,ier) d1 = dk1*tt*tr*d1zscale d2zscale = - chicdf(xzql,dfmp2,ier) + chicdf(xzqu,dfmp2,ier) d2 = dk2*tt*tr**2*d2zscale else d0 = 0.0d0 d1 = 0.0d0 d2 = 0.0d0 endif call matpv(wi,mmax,u(1,j),mmax,tempv,mmaxloc,m,m,1) do 620 i = 1,m temp(1,i) = tempv(i) 620 continue call matpv(tempv,mmaxloc,temp,1,tempm,mmaxloc,m,1,m) do 630 jj = 1,m do 631 k = 1,m tempm(k,jj) = (d2*tempm(k,jj)-d0*wi(k,jj)) / 2.0d0 631 continue tempv(jj) = d1*tempv(jj) 630 continue hu(1,1) = hu(1,1)+d0 do 640 jj=1,m tempvjj = tempv(jj) hu(1,jj+1) = hu(1,jj+1)+tempvjj hu(jj+1,1) = hu(jj+1,1)+tempvjj do 641 i=1,m hu(i+1,jj+1) = hu(i+1,jj+1)+tempm(i,jj) 641 continue 640 continue if ( (hsec1990()-timelim) .gt. hseclim ) then c give up ir = j go to 999 endif 600 continue 999 continue r=ir rscale= 1.0d0 / r do 700 j = 1,mp1 do 701 i = 1,mp1 hu(i,j) = hu(i,j)*rscale 701 continue 700 continue p(1) = hu(1,1) if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 800 j=1,mp1 do 801 i=1,mp1 hc(i,j)=hu(i,j)*pscale 801 continue 800 continue else do 900 j=1,mp1 do 901 i=1,mp1 hc(i,j)=-999.d0 901 continue 900 continue endif return end subroutine exaf(mmax,m,a,b,xmu,w,xl,p0,dm,dl) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c approximate answers by numerical quadrature. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension a(mmax),b(mmax),xmu(mmax),xl(mmax),w(mmax,mmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / s(mmaxloc),ss(mmaxloc,mmaxloc), / xmatmm(mmaxloc,mmaxloc),rmat(mmaxloc,mmaxloc), / acf(mmaxloc),bcf(mmaxloc), / xmup(mmaxloc),xlp(mmaxloc),wp(mmaxloc,mmaxloc) ccc / xlzt(1,mmaxloc), ccc / xlpzt(1,mmaxloc),xlmat(mmaxloc,mmaxloc), ccc / xlpmat(mmaxloc,mmaxloc) c mm1 = m - 1 do 100 i = 1,mm1 xmup(i) = xmu(i) xli = xl(i) xlp(i) = xli ccc xlzt(1,i) = xli ccc xlpzt(1,i) = xli 100 continue xmup(m) = xmu(m) + 1.d-3 xlm = xl(m) ccc xlzt(1,m) = xlm xlm = xlm + 1.d-3 xlp(m) = xlm ccc xlpzt(1,m) = xlm ccc call matpv(xl,mmax,xlzt,1,xlmat,mmaxloc,m,1,m) ccc call matpv(xlp,mmaxloc,xlpzt,1,xlpmat,mmaxloc,m,1,m) do 201 j = 1,m do 202 i = 1,m ccc wp(i,j) = w(i,j) - xlmat(i,j) + xlpmat(i,j) wp(i,j) = w(i,j) - xl(i)*xl(j) + xlp(i)*xlp(j) 202 continue 201 continue call pf(mmaxloc,m,a,b,xmu ,w ,p0) call pf(mmaxloc,m,a,b,xmup,w ,p1m) call pf(mmaxloc,m,a,b,xmu ,wp,p1w) call pf(mmaxloc,m,a,b,xmup,w ,dm) call pf(mmaxloc,m,a,b,xmu ,wp,dl) if(p0.ne.-999.d0.and.p1m.ne.-999.d0) then dm = (dm-p0)/1.d-3 else dm=-999.d0 endif if(p0.ne.-999.d0.and.p1w.ne.-999.d0) then dl = (dl-p0)/1.d-3 else dl=-999.d0 endif c extra test inserted by vah, sep 22, 92, to allow for inaccuracies c in the fortran pheev3 code (used inside pf): if(p0.eq.0.d0.and.p1m.eq.0.d0) then dm=-999.d0 endif if(p0.eq.0.d0.and.p1w.eq.0.d0) then dl=-999.d0 endif return end subroutine ffact(mmax,m,x,a,b,xmu,d,xl,y0,y1,y2) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c integrand in one-factor model. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) double precision nrmpdf dimension x(40),a(mmax),b(mmax),xmu(mmax),d(mmax),xl(mmax), / y0(40),y1(40),y2(40) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // xt(1,40),xxx(mmaxloc,40) c mm1 = m - 1 do 10 i = 1,40 xt(1,i) = x(i) 10 continue call matpv(xl,mmax,xt,1,xxx,mmax,m,1,40) do 101 j = 1,40 tt = 1.0d0 do 102 i = 1,mm1 xmui = xmu(i) xxxij = xxx(i,j) di = d(i) xhs = (b(i)-xmui-xxxij) / di rhs = (a(i)-xmui-xxxij) / di tt = tt*(pheev(xhs) - pheev(rhs)) 102 continue xmum = xmu(m) xxxmj = xxx(m,j) dm = d(m) xhs = (b(m)-xmum-xxxmj) / dm rhs = (a(m)-xmum-xxxmj) / dm y0(j) = tt*(pheev(xhs) - pheev(rhs)) ttt = nrmpdf(xhs,0.0d0,1.0d0) - nrmpdf(rhs,0.0d0,1.0d0) y1(j) = -tt*ttt / dm y2(j) = y1(j)*x(j) 101 continue return end subroutine fdl(mmax,m,h,xl,d) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c return specific elements from the H matrix. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension h(mmax+1,mmax+1),xl(mmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c mp1 = m + 1 d = 0.0d0 do 10 i = 1,m d = d + h(mp1,i+1)*xl(i) 10 continue d = d * 2.0d0 return end subroutine fghk(ii,mmax,m,xmu,w,wi,chw,a,b,ir,iseed,pm,p,hu) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c see GHK for documentation. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / p(1),hu(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / f(mmaxloc-1),g(mmaxloc-1,mmaxloc-1),vt(mmaxloc-1), / a2(mmaxloc-1),b2(mmaxloc-1), ccc--local variables for fghk, previous ones are from aus which calls it. / u(mmaxloc,irmaxloc), / tt(mmaxloc,irmaxloc),ttt(irmaxloc,mmaxloc), / wgt(1,irmaxloc),chw1(1,mmaxloc), / v(mmaxloc,irmaxloc),chw1tt(1,irmaxloc), / xmatmm(mmaxloc,mmaxloc),vv(mmaxloc) c mmaxp1 = mmax + 1 mp1 = m + 1 r = ir rscale = 1.0d0/r do 10 j = 1,ir do 11 i = 1,m cvahfix u(i,j) = unifrm(iseed) u(i,j) = randu(iseed) 11 continue 10 continue xmu1 = xmu(1) chw11 = chw(1,1) tascale = pheev((a(1)-xmu1)/(chw11+1.d-100)) tbscale = pheev((b(1)-xmu1)/(chw11+1.d-100)) do 20 j = 1,ir u1j = u(1,j) ccc should never happen: ccc if(u1j*tascale+(1.0d0-u1j)*tbscale.le.0.d0) then ccc print *,'error! phinv in FGHK (AUS/SUS)' ccc endif tt(1,j) = phinv(u1j*tascale+(1.0d0-u1j)*tbscale) pa = tbscale - tascale wgt(1,j) = pa 20 continue do 30 i = 2,m im1 = i - 1 do 31 k = 1,im1 chw1(1,k) = chw(i,k) 31 continue call matpv(chw1,1,tt,mmaxloc,chw1tt,1,1,im1,ir) xmui = xmu(i) chwii = chw(i,i) do 32 j = 1,ir chw1tt1j = chw1tt(1,j) uij = u(i,j) ta = pheev((a(i)-xmui-chw1tt1j)/(chwii+1.d-100)) tb = pheev((b(i)-xmui-chw1tt1j)/(chwii+1.d-100)) ccc should never happen: ccc if(uij*ta+(1.0d0-uij)*tb.le.0.d0) then ccc print *,'error! phinv in FGHK (AUS/SUS)' ccc endif tt(i,j) = phinv(uij*ta+(1.0d0-uij)*tb) wgt(1,j) = wgt(1,j)*(tb-ta) 32 continue 30 continue sum = 0.0d0 do 40 i = 1,ir sum = sum + wgt(1,i) 40 continue p(1) = sum/r if (ii .eq. 1) then do 45 j = 1,mp1 do 46 i = 1,mp1 hu(i,j) = 0.0d0 46 continue 45 continue else call matpv(chw,mmax,tt,mmaxloc,v,mmaxloc,m,m,ir) do 50 j=1,m sum = 0.0d0 do 51 i=1,ir val = v(j,i)*wgt(1,i) sum = sum + val ttt(i,j) = val 51 continue vv(j) = sum 50 continue call matpv(wi,mmax,v,mmaxloc,tt,mmaxloc,m,m,ir) call matpv(ttt,irmaxloc,wi,mmax,ttt,irmaxloc,ir,m,m) call matpv(tt,mmaxloc,ttt,irmaxloc,xmatmm,mmaxloc,m,ir,m) call matpv(wi,mmax,vv,mmaxloc,vv,mmaxloc,m,m,1) hu(1,1) = p(1) do 60 j=1,m vvj = vv(j) hu(1,j+1) = vvj*rscale hu(j+1,1) = vvj*rscale do 61 i=1,m hu(i+1,j+1) = (xmatmm(i,j)*rscale - p(1)*wi(i,j))*0.5d0 61 continue 60 continue endif pm = pa return end subroutine fi(mmax,irmax,m,ir,a,b,dkap,fi0,fi1,fi2) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c FI is a proc that returns partial moments of the standard normal. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension a(mmax,irmax),b(mmax,irmax),dkap(mmax,irmax), / fi0(mmax,irmax),fi1(mmax,irmax),fi2(mmax,irmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c do 101 j = 1,ir do 102 i = 1,m aij = a(i,j) bij = b(i,j) dkapij = dkap(i,j) fi0ij = pheev(bij) - pheev(aij) fi1ij = -dkapij*fi0ij - fi0ij fi2ij = fi0ij - dkapij*fi1ij - (bij - dkapij)*pheev(bij) * + (aij - dkapij)*pheev(aij) fi0(i,j) = fi0ij fi1(i,j) = fi1ij fi2(i,j) = fi2ij 102 continue 101 continue return end subroutine ghk(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c M. Geweke-Hajivassiliou-Keane Simulator (GHK) c The GHK procedure assumes that U is an array of uniform [0,1] c variates, independent within columns, and in general independent c between columns. It does not use PARM. It returns both HU and HC. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,m),wi(mmax,m),chw(mmax,m), / u(mmax,ir),parm(irmax+1,2), / p(1),hu(mmax+1,m+1),hc(mmax+1,m+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / tt(mmaxloc,irmaxloc),wgt(irmaxloc), / chw1(1,mmaxloc),chw1tt(1,irmaxloc), / ttt(irmaxloc,mmaxloc),v(mmaxloc,irmaxloc), / xmatmm(mmaxloc,mmaxloc),vdd(mmaxloc) c mp1 = m+1 mmaxp1 = mmax + 1 r = ir rscale = 1.0d0/r xmu1 = xmu(1) chw11 = chw(1,1) tascale = pheev((a(1)-xmu1)/(chw11+1.d-100)) tbscale = pheev((b(1)-xmu1)/(chw11+1.d-100)) do 10 j = 1,ir u1j = u(1,j) ccc should never happen: ccc if(u1j*tascale+(1.0d0-u1j)*tbscale.le.0.d0) then ccc print *,'GHK: warning! argument to phinv LE 0' ccc print *,u1j,tascale,tbscale ccc print *,u1j*tascale+(1.0d0-u1j)*tbscale ccc endif tt(1,j) = phinv(u1j*tascale+(1.0d0-u1j)*tbscale) wgt(j) = tbscale - tascale 10 continue do 20 i = 2, m im1 = i-1 do 21 k = 1,im1 chw1(1,k) = chw(i,k) 21 continue call matpv(chw1,1,tt,mmaxloc,chw1tt,1,1,im1,ir) xmui = xmu(i) chwii = chw(i,i) do 22 j = 1,ir chw1tt1j = chw1tt(1,j) uij = u(i,j) ta = pheev((a(i)-xmui-chw1tt1j)/(chwii+1.d-100)) tb = pheev((b(i)-xmui-chw1tt1j)/(chwii+1.d-100)) tt(i,j) = phinv(uij*ta+(1.0d0-uij)*tb) wgt(j) = wgt(j) * (tb-ta) 22 continue 20 continue call matpv(chw,mmax,tt,mmaxloc,v,mmaxloc,m,m,ir) sum = 0.0d0 do 25 j = 1,ir sum = sum + wgt(j) 25 continue p(1) = sum*rscale do 30 j = 1,m sumt = 0.0d0 do 31 i = 1,ir val = v(j,i)*wgt(i) sumt = sumt + val ttt(i,j) = val 31 continue vdd(j) = sumt 30 continue call matpv(wi,mmax,v,mmaxloc,tt,mmaxloc,m,m,ir) call matpv(ttt,irmaxloc,wi,mmax,ttt,irmaxloc,ir,m,m) call matpv(tt,mmaxloc,ttt,irmaxloc,xmatmm,mmaxloc,m,ir,m) call matpv(wi,mmax,vdd,mmaxloc,vdd,mmaxloc,m,m,1) hu(1,1) = p(1) do 40 j=1,m vddj = vdd(j) hu(1,j+1) = vddj*rscale hu(j+1,1) = vddj*rscale do 41 i=1,m hu(i+1,j+1) = (xmatmm(i,j)*rscale - p(1)*wi(i,j))*0.5d0 41 continue 40 continue if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 50 j=1,mp1 do 51 i=1,mp1 hc(i,j)=hu(i,j)*pscale 51 continue 50 continue else do 60 j=1,mp1 do 61 i=1,mp1 hc(i,j)=-999.d0 61 continue 60 continue endif return end subroutine gss(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c Q. Gibbs's Sampler Simulator (GSS) c The GSS simulates HC, but does not simulate P or HU. It assumes c that PARM[1,1] gives the number of rounds to be used in the c construction of the random draws (JL) and that U is an MxJL array of c uniform [0,1] random numbers. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(1), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter (mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / v(mmaxloc),wimkk(mmaxloc-1),wmkk(mmaxloc-1), / wmkmki(mmaxloc-1,mmaxloc-1),vbar(mmaxloc), / vvpbar(mmaxloc,mmaxloc),condsdk(mmaxloc), / condmukf(mmaxloc,mmaxloc-1),ui(mmaxloc,irmaxloc), / mkind(mmaxloc-1) c common/times/hseclim c mp1 = m + 1 mm1 = m - 1 r = ir jl = parm(1) do 10 j = 1,jl do 11 i=1,m ui(i,j)=u(i,j) 11 continue 10 continue do 100 ii = 1,m v(ii) = (b(ii)+a(ii))/2.0d0 - xmu(ii) 100 continue timestr = hsec1990() do 200 k=1,m if(k.eq.1) then do 201 n=1,mm1 mkind(n)=n+1 201 continue elseif(k.eq.m) then do 202 n=1,mm1 mkind(n)=n 202 continue else do 203 n=1,k-1 mkind(n)=n 203 continue do 204 n=k,mm1 mkind(n)=n+1 204 continue endif do 300 ni=1,mm1 wmkk(ni)=w(mkind(ni),k) wimkk(ni)=wi(mkind(ni),k) 300 continue do 301 nj=1,mm1 do 302 ni=1,mm1 wmkmki(ni,nj)=wi(mkind(ni),mkind(nj)) * - wimkk(ni)*wimkk(nj)/(wi(k,k)+1.d-100) 302 continue 301 continue temp=0.d0 do 303 nj=1,mm1 do 304 ni=1,mm1 temp=temp+wmkk(ni)*wmkmki(ni,nj)*wmkk(nj) 304 continue 303 continue condsdk(k)=dsqrt(w(k,k)-temp) do 305 nj=1,mm1 temp=0.d0 do 306 ni=1,mm1 temp=temp+wmkk(ni)*wmkmki(ni,nj) 306 continue condmukf(k,nj)=temp 305 continue 200 continue do 4001 nj=1,m vbar(nj) =0.d0 do 4002 ni=1,m vvpbar(ni,nj)=0.d0 4002 continue 4001 continue do 400 irepl=1,ir if (irepl.gt.1) then iseedgss=ui(1,1)*1.e+9 do 401 nj=1,jl do 402 ni=1,m ui(ni,nj)=randu(iseedgss) 402 continue 401 continue endif do 500 j=1,jl do 600 k=1,m if(k.eq.1) then do 2001 n=1,mm1 mkind(n)=n+1 2001 continue elseif(k.eq.m) then do 2002 n=1,mm1 mkind(n)=n 2002 continue else do 2003 n=1,k-1 mkind(n)=n 2003 continue do 2004 n=k,mm1 mkind(n)=n+1 2004 continue endif csdk=condsdk(k) cmuk=0.d0 do 403 nj=1,mm1 mk=mkind(nj) ccc print *,'nj,mk,k,mmaxloc',nj,mk,k,mmaxloc cerror??? cmuk=cmuk+condmukf(k,mk)*v(mk) cmuk=cmuk+condmukf(k,nj)*v(mk) 403 continue ta=pheev((a(k)-xmu(k)-cmuk)/(csdk+1.d-100)) tb=pheev((b(k)-xmu(k)-cmuk)/(csdk+1.d-100)) uikj=ui(k,j) v(k)=cmuk+csdk*phinv(uikj*tb+(1.d0-uikj)*ta+1.d-100) 600 continue 500 continue do 501 nj=1,m vbar(nj)=vbar(nj)+v(nj) do 502 ni=1,m vvpbar(ni,nj)=vvpbar(ni,nj)+v(ni)*v(nj) 502 continue 501 continue if((hsec1990()-timestr) .gt. hseclim) then c give up because of time-limit ir=i r=ir goto 4000 endif 400 continue 4000 continue c vah: potential to exploit SYMMETRY hc(1,1)=1.d0 do 700 nj=1,m temp=0.d0 do 701 i=1,m temp=temp+vbar(i)*wi(i,nj) 701 continue temp=temp/r hc(1,nj+1)=temp hc(nj+1,1)=temp do 702 ni=1,m res=0.d0 do 703 k=1,m do 704 j=1,m res=res+wi(ni,j)*(vvpbar(j,k)/r-w(j,k))*wi(k,nj)/2.d0 704 continue 703 continue hc(ni+1,nj+1)=res 702 continue 700 continue do 800 nj=1,mp1 do 801 ni=1,mp1 hu(ni,nj)=-999.d0 801 continue 800 continue p(1)=-999.d0 return end subroutine nden(mmax,irmax,m,ir,vd,wi,chw,q) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c NDEN is the multivariate normal density. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension vd(mmax,irmax),wi(mmax,mmax),chw(mmax,mmax),q(irmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // e(mmaxloc,irmaxloc) c dkk = 0.39894228040143 ** m prod = 1.0d0 do 10 i = 1,m prod = prod * chw(i,i) 10 continue dkk = dkk / prod call matpv(wi,mmax,vd,mmax,e,mmaxloc,m,m,ir) do 100 j = 1,ir sum = 0.0d0 do 101 i = 1,m sum = sum + vd(i,j)*e(i,j) 101 continue q(j)= dkk * dexp(-0.5d0*sum) 100 continue return end subroutine ndenist(mmax,irmax,m,ir,vd,wi,chw,q) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension vd(mmax,irmax),wi(mmax,mmax),chw(mmax,mmax),q(irmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / e(mmaxloc,irmaxloc),s(mmaxloc),vdc(mmaxloc,irmaxloc), / ymat(mmaxloc,irmaxloc) c dkk = 0.39894228040143 ** m prod = 1.0d0 do 10 i = 1,m prod = prod * chw(i,i) 10 continue dkk = dkk / (prod+1.d-100) call matpv(wi,mmax,vd,mmax,e,mmaxloc,m,m,ir) do 100 j = 1,ir sum = 0.0d0 do 101 i = 1,m sum = sum + vd(i,j)*e(i,j) 101 continue q(j)= dkk * dexp(-0.5d0*sum) 100 continue return end subroutine nise(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c J. Normal Importance Sampling Simulator (NIS) c This simulator provides HC = HU/P. It provides two alternative c simulators, corresponding to the exponential (NISE) and independent c truncated normal (NIST) comparison distributions. These simulators c assume U is an array of uniform random (0,1) variates. These procs c call the routine NDEN that returns a vector of values of the c multivariate normal density. c Procedure NISE does not use C or PARM. It calls the bilateral c exponential CDF BEXP and its inverse BEXPI from LIB_A.G. It c guarantees P positive. Provides both HU and HC. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,ir),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / e(mmaxloc,irmaxloc),s(mmaxloc), / diagw(mmaxloc),bxmu(mmaxloc),axmu(mmaxloc), / t1(mmaxloc),t2(mmaxloc),t3(mmaxloc),t4(mmaxloc), / temp(mmaxloc),pa(mmaxloc,irmaxloc),pb(mmaxloc,irmaxloc), / wgt(irmaxloc),xmatirm(irmaxloc,mmaxloc),q(irmaxloc), / humm(mmaxloc,mmaxloc),vdd(mmaxloc),vd(mmaxloc,irmaxloc), / d(mmaxloc,irmaxloc) c mp1 = m + 1 do 10 i = 1,m diagw(i) = w(i,i) s(i) = dsqrt(diagw(i)) bxmu(i) = b(i) - xmu(i) axmu(i) = a(i) - xmu(i) temp(i) = -3.0d0 * s(i) q(i) = 3.0d0 * s(i) vdd(i) = 0.0d0 10 continue do 20 i = 1,m t1(i) = 0.0d0 t2(i) = 0.0d0 t3(i) = 0.0d0 t4(i) = 0.0d0 if (bxmu(i).lt.temp(i)) t1(i)=1.0d0 if (axmu(i).gt.q(i)) t2(i)=1.0d0 if (bxmu(i).lt.q(i)) t3(i)=1.0d0 if (axmu(i).gt.temp(i)) t4(i)=1.0d0 t1(i) = t1(i)*(bxmu(i)-2.0d0*s(i)) * + t2(i)*axmu(i) * + (1.0d0-t1(i)-t2(i))* * (temp(i)+t4(i)*(axmu(i)+q(i))) t2(i) = t1(i)*bxmu(i) * + t2(i)*(axmu(i)+2.0d0*s(i)) * + (1.0d0-t1(i)-t2(i))* * (q(i)+t3(i)*(bxmu(i)+temp(i))) 20 continue do 30 i = 1,m t3(i) = 0.0d0 t4(i) = 0.0d0 if (t2(i).lt.vdd(i)) t4(i)=1.0d0 if (t1(i).gt.vdd(i)) t3(i)=1.0d0 s(i) = 1.0d0 - t4(i) - t3(i) temp(i) = 1.d-100 * + (t2(i)-t1(i))*(0.5d0-s(i)/4.0d0)/(diagw(i)+1.d-100) 30 continue call bexp(mmaxloc,m,axmu,temp,axmu) call bexp(mmaxloc,m,bxmu,temp,bxmu) do 40 j = 1,ir do 41 i =1,m pa(i,j) = axmu(i) pb(i,j) = bxmu(i) d(i,j) = u(i,j)*pa(i,j) + (1.0d0-u(i,j))*pb(i,j) 41 continue 40 continue r = ir rscale = 1.0d0/r call bexpi(mmaxloc,m,ir,d,temp,pa,pb,vd,d) call nden(mmaxloc,irmaxloc,m,ir,vd,wi,chw,q) sum = 0.0d0 do 50 j = 1,ir prod = 1.0d0 do 51 i = 1,m prod = prod*d(i,j) 51 continue wgt(j) = q(j) / (prod + 1.0d-100) sum = sum + wgt(j) do 52 i = 1,m pa(i,j) = wgt(j)*vd(i,j) 52 continue 50 continue p(1) = sum*rscale call matpv(wi,mmax,pa,mmaxloc,pb,mmaxloc,m,m,ir) call matptv(vd,mmaxloc,wi,mmax,xmatirm,irmaxloc,m,ir,m) call matpv(pb,mmaxloc,xmatirm,irmaxloc,humm,mmaxloc,m,ir,m) do 60 j = 1, m do 61 i = 1,m humm(i,j) = humm(i,j)*rscale - p(1)*w(i,j) 61 continue 60 continue do 70 j = 1,m sum = 0.0d0 do 71 i = 1,ir sum = sum + pa(j,i) 71 continue vdd(j) = sum*rscale 70 continue call matpv(wi,mmax,vdd,mmaxloc,vdd,mmaxloc,m,m,1) hu(1,1) = p(1) do 80 j=1,m hu(1,j+1) = vdd(j) hu(j+1,1) = vdd(j) do 81 i=1,m hu(i+1,j+1) = humm(i,j)*0.5d0 81 continue 80 continue if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 90 j=1,mp1 do 91 i=1,mp1 hc(i,j)=hu(i,j)*pscale 91 continue 90 continue else do 100 j=1,mp1 do 101 i=1,mp1 hc(i,j)=-999.d0 101 continue 100 continue endif return end subroutine nist(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c see NISE for documentation. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / e(mmaxloc,irmaxloc),s(mmaxloc),vdc(mmaxloc,irmaxloc), / ymat(mmaxloc,irmaxloc),xmatirm(irmaxloc,mmaxloc), / y(irmaxloc),wgt(irmaxloc),vd(mmaxloc,irmaxloc), / temp(mmaxloc),humm(mmaxloc,mmaxloc),vdd(mmaxloc) c mp1 = m + 1 r = ir rscale = 1.0d0/r do 10 i = 1,m temp(i) = a(i) - xmu(i) vdd(i) = b(i) - xmu(i) 10 continue call rndt(mmax,irmax,m,ir,u,w,temp,vdd,vd,y) call ndenist(mmax,irmax,m,ir,vd,wi,chw,wgt) sum = 0.0d0 do 20 j = 1,ir wgt(j) = wgt(j)/(y(j)+1.0d-100) sum = sum + wgt(j) do 21 i = 1,m ymat(i,j) = wgt(j)*vd(i,j) 21 continue 20 continue p(1) = sum*rscale do 30 j = 1,m sum = 0.0d0 do 31 i = 1,ir sum = sum + ymat(j,i) 31 continue vdd(j) = sum*rscale 30 continue call matpv(wi,mmax,ymat,mmaxloc,ymat,mmaxloc,m,m,ir) call matptv(vd,mmaxloc,wi,mmax,xmatirm,irmaxloc,m,ir,m) call matpv(ymat,mmaxloc,xmatirm,irmaxloc,humm,mmaxloc,m,ir,m) do 40 j=1,m do 41 i=1,m humm(i,j) = humm(i,j)*rscale - wi(i,j)*p(1) 41 continue 40 continue call matpv(wi,mmax,vdd,mmaxloc,vdd,mmaxloc,m,m,1) hu(1,1) = p(1) do 50 j=1,m hu(1,j+1) = vdd(j) hu(j+1,1) = vdd(j) do 51 i=1,m hu(i+1,j+1) = humm(i,j)*0.5d0 51 continue 50 continue if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 60 j=1,mp1 do 61 i=1,mp1 hc(i,j)=hu(i,j)*pscale 61 continue 60 continue else do 70 j=1,mp1 do 71 i=1,mp1 hc(i,j)=-999.d0 71 continue 70 continue endif return end subroutine pcf(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c N. Parabolic Cylinder Function Simulator (PCF) c The PCF proc assumes that U is an array whose column vectors are c points in the intersection of the unit sphere and the negative c orthant. c The routine assumes that PARM(1,1) = DET(W). c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / aa(mmaxloc,irmaxloc),bb(mmaxloc), / vwv(irmaxloc),bwv(irmaxloc), / svec(irmaxloc),ro(irmaxloc), / y0(irmaxloc),y1(irmaxloc),y2(irmaxloc),tempv(mmaxloc) c mp1 = m + 1 r = ir c cvah det(w) must be passed as parm(1,1) detw=parm(1,1) c c dkp dm = m dkp = dgamma(dm/2.0d0) * (2.0d0**(dm/2.0d0-1.0d0)) * dsqrt(detw) dkp = 1.0d0/(dkp + 1.d-100) c c bb=wi*(mu-b) (mx1) do 20 i=1,m sum=0.d0 do 21 k=1,m sum = sum + wi(i,k)*(xmu(k) - b(k)) 21 continue bb(i) = sum 20 continue c c aa=wi*u (mxr) call matpv(wi,mmax,u,mmax,aa,mmaxloc,m,m,ir) c c vwv=sumc(u.*aa) (rx1) do 30 j = 1,ir sum = 0.0d0 do 31 k = 1,m sum = sum + u(k,j)*aa(k,j) 31 continue vwv(j) = sum 30 continue c c bwv=u'bb (rx1) do 40 j=1,ir sum=0.d0 do 41 k=1,m sum = sum + u(k,j)*bb(k) 41 continue bwv(j) = sum 40 continue c c svec, ro c sval=-(xmu-b)'bb/2 sval=0.d0 do 50 k=1,m sval = sval - (xmu(k)-b(k))*bb(k) 50 continue sval = sval/2.d0 do 51 j=1,ir c svec temp = sval+bwv(j)**2/(2.d0*vwv(j)+1.d-100) svec(j) = dkp*dexp(temp) c ro valmin = (a(1)-b(1))/(u(1,j)+1.d-100) do 52 i=2,m temp = (a(i)-b(i))/(u(i,j)+1.d-100) if (temp.lt.valmin) valmin = temp 52 continue ro(j) = valmin 51 continue c c y0, y1, y2 c note slight difference from gauss: k=m-1 is calculated INSIDE cf: call cf(irmaxloc,m,ir,vwv,bwv,ro,y0,y1,y2) c p sum = 0.0d0 do 55 i = 1,ir sum = sum + y0(i)*svec(i) 55 continue p(1) = sum/r c c tempv do 60 i = 1,m sum = 0.d0 do 61 j = 1,m do 62 k = 1,ir sum = sum + wi(i,j)*u(j,k)*y1(k)*svec(k) 62 continue 61 continue tempv(i) = sum/r 60 continue c c hu twotom=2.d0**m twotomp1=twotom*2.d0 hu(1,1) = p(1)/twotom do 70 j=1,m temp = (tempv(j) - p(1)*bb(j)) / twotom hu(1,j+1) = temp hu(j+1,1) = temp c should exploit SYMMETRY!!! do 71 i=1,m sum = 0.d0 do 72 k=1,ir sum = sum + aa(i,k)*aa(k,i)*y2(k)*svec(k) 72 continue temp = * (- tempv(i)*bb(j) - tempv(j)*bb(i) * + sum/r + p(1)*(bb(i)*bb(j) - wi(i,j)) ) / twotomp1 hu(i+1,j+1) = temp 71 continue 70 continue c c hc p(1) = hu(1,1) if ( p(1) .gt. 0.0d0 ) then do 80 j = 1,mp1 do 81 i =1,mp1 hc(i,j) = hu(i,j)/p(1) 81 continue 80 continue else do 90 j = 1,mp1 do 91 i = 1,mp1 hc(i,j) = -999.0d0 91 continue 90 continue endif return end subroutine pf(mmax,m,aa,bb,xmu,w,p) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c used by EXAF. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension aa(mmax),bb(mmax),xmu(mmax),w(mmax,mmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / s(mmaxloc),ss(mmaxloc,mmaxloc), / xmatmm(mmaxloc,mmaxloc),rmat(mmaxloc,mmaxloc), / a(mmaxloc),b(mmaxloc) c do 100 i = 1,m s(i) = dsqrt(w(i,i)) ss(i,i) = 1.0d0/s(i) 100 continue call matpv(ss,mmaxloc,w,mmax,xmatmm,mmaxloc,m,m,m) call matpv(xmatmm,mmaxloc,ss,mmaxloc,rmat,mmaxloc,m,m,m) do 110 i = 1,m a(i) = (aa(i)-xmu(i)) / s(i) b(i) = (bb(i)-xmu(i)) / s(i) 110 continue if (m.eq.2) then r = rmat(1,2) a1 = a(1) a2 = a(2) b1 = b(1) b2 = b(2) p = pheev2(b1,b2,r) - pheev2(a1,b2,r) * - pheev2(b1,a2,r) + pheev2(a1,a2,r) else if (m.eq.3) then r1 = rmat(1,2) r2 = rmat(2,3) r3 = rmat(3,1) a1 = a(1) a2 = a(2) a3 = a(3) b1 = b(1) b2 = b(2) b3 = b(3) p = pheev3(b1,b2,b3,r1,r2,r3) - pheev3(a1,b2,b3,r1,r2,r3) p = p - pheev3(b1,a2,b3,r1,r2,r3) - pheev3(b1,b2,a3,r1,r2,r3) p = p + pheev3(a1,a2,b3,r1,r2,r3) + pheev3(b1,a2,a3,r1,r2,r3) p = p + pheev3(a1,b2,a3,r1,r2,r3) - pheev3(a1,a2,a3,r1,r2,r3) else p = -999.0d0 endif return end subroutine pker(z,mmax,irmax,m,ir,prodzwgt) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c see PKF for documentation. c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Polynomial kernel function; PKFS calls this routine. c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension z(mmax,irmax),prodzwgt(irmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c c common // c / onesmir(mmaxloc,irmaxloc), c / zerosmir(mmaxloc,irmaxloc),xm1mir(mmaxloc,irmaxloc) c c do 10 j = 1,ir c do 11 i = 1,m c onesmir(i,j) = 1.0d0 c zerosmir(i,j) = 0.0d0 c xm1mir(i,j) = - 1.0d0 c 11 continue c 10 continue c call hadcomp(z,mmax,onesmir,mmaxloc,onesmir,mmaxloc,m,ir,'lt') c call hadcomp(z,mmax,zerosmir,mmaxloc,zerosmir,mmaxloc,m,ir,'lt') c call hadcomp(z,mmax,xm1mir,mmaxloc,xm1mir,mmaxloc,m,ir,'lt') do 101 j = 1,ir prod = 1.0d0 do 102 i = 1,m z1 = 0.0d0 z0 = 0.0d0 zm1 = 0.0d0 if (z(i,j).lt.(1.0d0)) z1=1.0d0 if (z(i,j).lt.(0.0d0)) z0=1.0d0 if (z(i,j).lt.(-1.0d0)) zm1=1.0d0 factor = 0.5d0*(1.0d0-z(i,j)*(2.0d0-z(i,j)))*z1 * - z(i,j)**2*z0 * + 0.5d0*( 1.0d0 + z(i,j)*(2.0d0+z(i,j)) )*zm1 * prod = prod*factor 102 continue prodzwgt(j) = prod 101 continue return end subroutine pkfs(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c K. Polynomial Kernel-Smoothed Frequency Simulator (PKF) c The KFS simulator uses the kernel (31), defined by the function c PKER. It assumes that U is an array of standard normal variates, c independent within rows. It does not use W. PARM[1,1] contains a c window width parameter. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / vd(mmaxloc,irmaxloc), / vdbb(mmaxloc,irmaxloc),vdaa(mmaxloc,irmaxloc), / pkerbb(irmaxloc),pkeraa(irmaxloc), / wivd(mmaxloc,irmaxloc),aat(irmaxloc,mmaxloc), / humm(mmaxloc,mmaxloc), / vdd(mmaxloc) c mmaxlcp1 = mmaxloc + 1 mp1 = m + 1 mmaxp1 = mmax + 1 r = ir rscale = 1.d0/r c call matpv(chw,mmax,u,mmax,vd,mmaxloc,m,m,ir) do 10 j = 1,ir do 11 i = 1,m vdbb(i,j) = (vd(i,j)-(b(i)-xmu(i)))/parm(1,1) vdaa(i,j) = (vd(i,j)-(a(i)-xmu(i)))/parm(1,1) 11 continue 10 continue call pker(vdbb,mmaxloc,irmaxloc,m,ir,pkerbb) call pker(vdaa,mmaxloc,irmaxloc,m,ir,pkeraa) call matpv(wi,mmax,vd,mmax,wivd,mmaxloc,m,m,ir) sum = 0.0d0 do 20 j = 1,ir fval = pkerbb(j) - pkeraa(j) sum = sum + fval do 21 i = 1,m aat(j,i) = wivd(i,j)*fval 21 continue 20 continue p(1) = sum*rscale call matpv(wivd,mmaxloc,aat,irmaxloc,humm,mmaxloc,m,ir,m) do 30 j=1,m sum = 0.0d0 do 31 i=1,ir sum = sum + aat(i,j) 31 continue vdd(j) = sum*rscale 30 continue hu(1,1) = p(1) do 40 j=1,m hu(1,j+1) = vdd(j) hu(j+1,1) = vdd(j) do 41 i=1,m hu(i+1,j+1) = (humm(i,j)*rscale - wi(i,j)*p(1))*0.5d0 41 continue 40 continue if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 50 j = 1,mp1 do 51 i = 1,mp1 hc(i,j) = hu(i,j)*pscale 51 continue 50 continue else do 60 j = 1,mp1 do 61 i = 1,mp1 hc(i,j) = -999.0d0 61 continue 60 continue endif return end subroutine reswrite(cfil,line,irec) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c direct access version c Reswrite writes the character string line to the end of the cfil. c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ character*(*) cfil,line open(10,file=cfil,access='direct',form='formatted',recl=130) ccc print *,'inside new: line',line write(10,'(a)',rec=irec) line(1:ltrim(line)) close(10) end subroutine rndbase(kmax,k,iseed1,v,t) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c RNDBASE is a routine that draws a random basis in K space, where K is c the dimension of the random coefficient vector. c RNDBASE takes as input a dimension K and returns a K by K array V that c is column orthonormal and random, bordered by a column of random c lengths of K-dim. normal random vectors. SEED1 is a seed for the c random number generator. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension v(kmax,k),t(kmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / vt(mmaxloc,mmaxloc),vvt(mmaxloc,mmaxloc) c do 100 j=1,k do 101 i=1,k v(i,j) = ranorm(0.0d0,1.0d0,iseed) vt(j,i) = v(i,j) 101 continue 100 continue call matpv(v,kmax,vt,mmaxloc,vvt,mmaxloc,k,k,k) 888 continue s = 0.0d0 do 110 i=1,k s = s + v(i,1)**2 110 continue if (s .lt. 1.d-16) then do 120 i=1,k v(i,1) = ranorm(0.0d0,1.0d0,iseed1) 120 continue go to 888 endif scale = 1.0d0 / (dsqrt(s)+1.d-100) do 200 i = 1,k v(i,1) = scale * v(i,1) t(i) = vvt(i,i) 200 continue do 300 kk = 2,k 999 continue do 310 i = 1,kk rhs = 0.0d0 do 311 j = 1,k rhs = rhs + v(j,kk)*v(j,i) v(j,kk) = v(j,kk) - v(j,i)*rhs 311 continue 310 continue s = 0.0d0 do 320 i = 1,k s = s + v(i,kk)**2 320 continue if (s.lt. 1.d-16) then call mrnorm(v(1,kk),kmax,k,1,iseed1) go to 999 endif scale = 1.0d0 / (dsqrt(s)+1.d-100) do 330 i = 1,k v(i,kk) = scale * v(i,kk) 330 continue 300 continue return end subroutine rndt(mmax,irmax,m,ir,u,w,a,b,v,y) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c RNDT is a proc for drawing an array from truncated independent c normals that have the same means and variances as the multivariate c density. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension a(mmax),b(mmax),w(mmax,mmax),u(mmax,irmax), / v(mmax,irmax),y(irmax) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / e(mmaxloc,irmaxloc),s(mmaxloc),vdc(mmaxloc,irmaxloc), / ymat(mmaxloc,irmaxloc) c do 10 j = 1,m do 10 i = 1,m if (i.eq.j) then s(i) = dsqrt(w(i,j)) + 1.d-100 endif 10 continue dkk = 0.39894228040143 do 40 j = 1,ir prod = 1.0d0 do 50 i = 1,m bs= pheev( b(i)/s(i) ) as= pheev( a(i)/s(i) ) vdc(i,j) = phinv( u(i,j)*bs+(1.0d0-u(i,j))*as ) v(i,j) = s(i)*vdc(i,j) ccc---truncated normal density sdk = dkk/s(i) ymat(i,j) = sdk*dexp(-0.5d0*vdc(i,j)*vdc(i,j)) * /(bs-as+1.d-100) prod = prod * ymat(i,j) 50 continue y(j) = prod 40 continue return end subroutine sds(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c L. Stern Decomposition Simulator (SDS) c The SDS proc assumes that PARM[1,1] contains a scalar l > 0 such c that W-lI is positive definite, and assumes that C is a Cholesky c factor of W-lI. It assumes that U is an array of standard normal c variates, independent within each column. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,irmax),parm(irmax+1,2), / p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c include 'hdim.inc' parameter(mmaxloc=mmaxim,irmaxloc=irmaxim) c common // / v(mmaxloc,irmaxloc),ta(mmaxloc,irmaxloc), / tb(mmaxloc,irmaxloc),tk(mmaxloc,irmaxloc), / fb0(mmaxloc,irmaxloc),fb1(mmaxloc,irmaxloc), / fb2(mmaxloc,irmaxloc),qi1(mmaxloc,irmaxloc), / qi2(mmaxloc,irmaxloc),s(mmaxloc,irmaxloc), / dd(mmaxloc,mmaxloc),tt(mmaxloc), / tempmm(mmaxloc,mmaxloc),ttt(1,mmaxloc) c mp1 = m + 1 mmaxp1 = mmax + 1 mmaxlcp1 = mmaxloc + 1 r = ir rscale = 1.0d0 / r xlambda = dsqrt(parm(1,1)) call matpv(chw,mmax,u,mmax,v,mmaxloc,m,m,ir) do 101 j = 1,ir do 102 i = 1,m tb(i,j) = (b(i)-xmu(i)-v(i,j)) / xlambda ta(i,j) = (a(i)-xmu(i)-v(i,j)) / xlambda tk(i,j) = -2.0d0*xmu(i) / xlambda 102 continue 101 continue call fi(mmaxloc,irmaxloc,m,ir,ta,tb,tk,fb0,fb1,fb2) do 110 j = 1,ir do 120 i = 1,m if ( fb0(i,j) .eq. 0.0d0 ) then fb0(i,j) = .000000001 endif 120 continue 110 continue do 201 j = 1,ir do 202 i = 1,m qi1(i,j) = xlambda*fb1(i,j) qi2(i,j) = xlambda**2*fb2(i,j) s(i,j) = qi1(i,j)/fb0(i,j) 202 continue 201 continue do 210 j = 1,mp1 do 211 i = 1,mp1 hu(i,j) = 0.0d0 211 continue 210 continue do 300 j = 1,ir qq = 1.0d0 do 310 jk=1,m qq = qq*fb0(jk,j) do 311 ik=1,m if (ik.eq.jk) then dd(ik,jk) = qi2(ik,j)/fb0(ik,j) else dd(ik,jk) = 0.0d0 endif 311 continue 310 continue call matpv(wi,mmax,s(1,j),mmaxloc,tt,mmaxloc,m,m,1) call matpv(wi,mmax,dd,mmax,tempmm,mmax,m,m,m) call matpv(tempmm,mmax,wi,mmax,tempmm,mmax,m,m,m) do 320 ik = 1,m ttt(1,ik) = tt(ik) 320 continue call matpv(tt,mmaxloc,ttt,1,dd,mmaxloc,m,1,m) hu(1,1) = hu(1,1)+qq scale = qq*0.5d0 do 330 jk=1,m hu(1,jk+1) = hu(1,jk+1)+tt(jk)*qq hu(jk+1,1) = hu(jk+1,1)+tt(jk)*qq do 331 i=1,m hu(i+1,jk+1) = hu(i+1,jk+1) * +(tempmm(i,jk)-wi(i,jk)+dd(i,jk))*scale 331 continue 330 continue 300 continue do 400 j = 1,mp1 do 401 i = 1,mp1 hu(i,j) = hu(i,j)*rscale 401 continue 400 continue p(1) = hu(1,1) if (p(1) .gt. 0.0d0) then pscale = 1.0d0/p(1) do 500 j=1,mp1 do 501 i=1,mp1 hc(i,j)=hu(i,j)*pscale 501 continue 500 continue else do 600 j=1,mp1 do 601 i=1,mp1 hc(i,j)=-999.d0 601 continue 600 continue endif return end subroutine sim(mmax,m,xl,e1,p,hu,hc,ec,pp,dm,xlm,dl,xll) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c procedure to accumulate results. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit double precision (a-h,o-z) implicit integer*4 (i-n) dimension xl(mmax),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c pp = p dm = hu(m+1,1) xlm = hc(m+1,1) if ( hu(m+1,1) .ne. -999.0d0 ) then call fdl(mmax,m,hu,xl,dl) else dl = -999.0d0 endif if ( hc(m+1,1) .ne. -999.0d0 ) then call fdl(mmax,m,hc,xl,xll) else xll = -999.0d0 endif e2 = hsec1990() ec = e2 - e1 e1 = e2 return end subroutine sus(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) c===============================================================c c written by Vassilis Hajivassiliou c c with the aid of Yoosoon Chang. c c version 1.0, March 22, 1992 c c version 1.1, July 11, 1992 c c===============================================================c c c R. Sequentially Unbiased Simulator (SUS) c The SUS simulator is obtained as the product of an unbiased c simulator of H and an unbiased simulator of 1/P. c Computationally, if PARM(2,2) is made large (say 1.e50), SUS is the c same as the Asymptotically Unbiased Simulator (AUS). c c NB: Upon exit, the first element of the BLANK common block, contains c -- the maximum number of CFS draws allowed in AUS c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c implicit double precision (a-h, o-z) implicit integer*4 (i-n) dimension xmu(mmax),a(mmax),b(mmax), / w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax), / u(mmax,ir),parm(irmax+1,2), / hu(mmax+1,mmax+1),hc(mmax+1,mmax+1) c c simply calls aus with a very large PARM(2,2) ccc parm(2,2) = 1.0d50 parm(2,2) = 1.0d7 call aus(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc) return end