c c--------------------------------------------------------------------- c c 9/29/1997 "m-simul.f" S. Stern c c--------------------------------------------------------------------- c c this file contains subroutines to simulate the distribution of c missing values of variables for an observation conditional on c the observed variables for that observation. it takes as input c the joint distribution of the variables. c c--------------------------------------------------------------------- c subroutine choles(amat,nmats,nmat,cmat) c c this subroutine finds the cholesky decomposition of amat c and stores it in cmat c implicit real*8(a-h,o-z) dimension amat(nmats),cmat(nmats),c(122,122) do 5 j=1,nmat do 5 i=1,j 5 c(i,j)=0. l=1 do 10 j=1,nmat do 10 i=1,j if(i.eq.j)then x=0. do 20 k=1,i-1 20 x=x+c(k,i)**2 x=amat(l)-x if(x.gt.0)c(i,i)=x**.5 endif if(i.ne.j)then if(c(i,i).gt.0.)then x=0. do 30 k=1,i-1 30 x=x+c(k,i)*c(k,j) c(i,j)=(amat(l)-x)/c(i,i) endif c if(c(i,i).le.0.)c(i,j)=0. if(c(i,i).lt.0.)then write(6,100) 100 format(1x,'bad covariance matrix') stop endif endif l=l+1 10 continue l=1 do 40 j=1,nmat do 40 i=1,j cmat(l)=c(i,j) 40 l=l+1 return end c c--------------------------------------------------------------------- c subroutine errget(e,nch,omegas,nchs,u) c c this subroutine generates errors distributed normally with c mean zero and covariance matrix omega (omegas is the cholesky c decomposition of omega) c implicit real*8(a-h,o-z) dimension omegas(nchs),e(nch),u(nch) c c transpose by omegas c do 2 i=1,nch e(i)=0. do 2 j=1,i ij=(((i-1)*i)/2)+j 2 e(i)=e(i)+(u(j)*omegas(ij)) return end c c--------------------------------------------------------------------- c subroutine nrandom(iseed,z1) implicit real(a-h,o-z) x1=rand2(iseed) if(x1.gt..9999999)x1=.9999999 if(x1.lt..0000001)x1=.0000001 call mdnris(x1,z1,ier) if(abs(z1).gt.19.999)then if(z1.gt.19.999)z1=19.999 if(z1.lt.-19.999)z1=-19.999 endif return end c c--------------------------------------------------------------------- c subroutine simul c c this is the driver. c parameter(nobss=1000,nxs=122,ncs=8,nx2s=nxs*2,nxs2=(nxs+1)*nxs/2) implicit real*8(a-h,o-z) character*8 aids(3) real arg,farg,ut dimension acovt(nxs,nxs),amat(nxs2),ameant(nxs),amnsim(nxs), & amntemp(nxs),anumocc(7),anumwks(5),asigsim(nxs,nxs),cmat(nxs2), & concov(nxs,nxs),conmean(nxs),covinv(nxs,nxs),covsym(nxs2), & covtemp(nxs,nxs),cross(nxs,nxs),crotemp(nxs),dev(nxs),devt(nxs), & e(nxs),eigval(nxs),eigvec(nxs,nxs),et(nxs),ialmis(nxs), & imcou(nxs),imcout(nxs,7), & itr(nxs),itrco(3),itrinv(nxs),jco(5),jtrans(5), & nmcou(7),nmt(nobss),nvehic(4),nxbnryt(nxs),nxorddt(nxs), & scov(nxs,nxs,3),smean(nxs,2),sob(nxs,2),sob2(nxs,nxs),std(3), & u(nxs),ucov(nxs,nxs),umean(nxs),work1(nxs),work2(nx2s), & x(nobss,nxs),xrealt(nxs),xt(nxs),xtemp(nobss),xtr(nxs), & ytemp(nobss) common/nshare/nid(nobss,3),imx(nobss,nxs),nxbnry(nobss,nxs), & nxordd(nobss,nxs),ids(nxs),nchat(nxs),iestflg,iacrj,itype, & nvar,nfac common/xshare/xreal(nobss,nxs),amean(nxs),acov(nxs,nxs), & chat(nxs,ncs) data aids/'binary','ordered','contin'/ data anumocc/0.,1.,2.,4.5,11.5,24.5,50.0/ data anumwks/0.,0.5,2.0,26.0,104.0/ data jtrans/4,5,1,2,3/ data nvehic/0,1,2,4/ c c--------------------------------------------------------------------- c c initialize. c isetm=0 c c run program on monte carlo data, if "imonte = 1". c cmonte imonte = 0 imonte = 1 if (imonte .eq. 0) then nobs=584 nx=122 endif if (imonte .eq. 1) then nobs=1000 charry nobs=250 nx=15 endif c do 77 i=1,nx if(ids(i).eq.1)then do 78 j=1,nobs 78 x(j,i)=nxbnry(j,i) endif if(ids(i).eq.2)then do 79 j=1,nobs 79 x(j,i)=nxordd(j,i) endif if(ids(i).eq.3)then do 80 j=1,nobs 80 x(j,i)=xreal(j,i) endif 77 continue d1=-1. nx2=(nx+1)*nx/2 ndraw=10 iopout=3 iouttype=2 istan=0 c c--------------------------------------------------------------------- c c do for each facility type: c do 52 ii=itype,itype c c--------------------------------------------------------------------- c c check for variables always missing. c do 91 i=1,nx ialmis(i)=1 do 92 j=1,nobs if(nid(j,1).eq.ii)then if(imx(j,i).eq.1)then ialmis(i)=0 goto 91 endif endif 92 continue 91 continue c c------------------------------------------------------------------------- c c read parameters. c if(iopout.eq.2)open(unit=18,file='draws.dat') if(iopout.eq.3)open(unit=18,file='m-simul.dat') c open(unit=12,file="xmean.dat") open(unit=13,file="xcut.dat") open(unit=14,file="xcov.dat") open(unit=15,file="m-simul.out") c read(12,209) (amean(k),k=1,nx) read(13,209) ((chat(k,l),l=1,ncs),k=1,nx) read(14,209) ((acov(k,l),l=1,nx),k=1,nx) 209 format(/,5(2x,g13.6)) kl=0 do 51 k=1,nx do 51 l=1,k kl=kl+1 51 covsym(kl)=acov(k,l) call eigrs(covsym,nx,0,eigval,eigvec,nxs,work1,ier) if(eigval(1).le.0.)then write(6,202)(eigval(k),k=1,nx) 202 format(1x,'negative eigenvalue:',/,122(1x,g15.8)) stop endif close(12) close(13) close(14) c write(15,220) 220 format(//,4x,'m-simul.f -- simulation begins',//) c c------------------------------------------------------------------------- c imonteg=0 if(imonteg.eq.1)then itrc=0 do 75 i=1,nx if(ids(i).eq.3)then itrc=itrc+1 itr(itrc)=i itrinv(i)=itrc endif 75 continue itrco(1)=itrc do 74 i=1,nx if(ids(i).eq.1)then itrc=itrc+1 itr(itrc)=i itrinv(i)=itrc endif 74 continue itrco(2)=itrc-itrco(1) do 63 i=1,nx if(ids(i).eq.2)then itrc=itrc+1 itr(itrc)=i itrinv(i)=itrc endif 63 continue itrco(3)=itrc-(itrco(1)+itrco(2)) write(6,205)(i,itr(i),i=1,nx) 205 format(1x,'itr: ',(1x,i3,2x,i3),121(/,7x,i3,2x,i3)) write(6,204)(i,itrinv(i),i=1,nx) 204 format(1x,'itrinv: ',(1x,i3,2x,i3),121(/,10x,i3,2x,i3)) open(unit=18,file='monte.p') write(18,*)(itrco(i),i=1,3),nobs write(18,208)(amean(itr(j)),j=1,nx) 208 format(25(1x,5(f8.4,2x),/)) do 62 i=1,nx itri=itr(i) 62 write(18,207)(acov(itri,itr(j)),j=1,i) 207 format(1x,5(f8.4,2x),/,30(11x,4(f8.4,2x),/)) if(itrco(1).gt.0)write(18,*)(1,j=1,itrco(1)) njunk1=itrco(1)+itrco(2)+1 if(itrco(3).gt.0)then write(18,206)((nchat(itr(j))-1),j=njunk1,nx) 206 format(5(1x,25(i1,1x),/)) do 61 i=njunk1,nx itri=itr(i) nchatrt=nchat(itri) 61 write(18,208)(chat(itri,j),j=3,nchatrt) endif write(18,208)(.15,i=1,nx) stop endif c c------------------------------------------------------------------------- c do 38 i=1,nx sob(i,1)=0. sob(i,2)=0. do 50 j=1,2 50 smean(i,j)=0. do 39 j=1,nx sob2(i,j)=0. do 39 k=1,3 39 scov(i,j,k)=0. 38 continue c c--------------------------------------------------------------------- c c do for each observation. c do 2 i=1,nobs c c check if facility is the correct type. c if(nid(i,1).eq.ii)then c c check if facility has all variables missing. c initialize type counters. c do 3 j=1,7 3 nmcou(j)=0 c c determine status of each variable. c do 4 j=1,nx if((imx(i,j).eq.0).and.(ialmis(j).eq.0))then if(ids(j).eq.1)then c c status is missing binary. c nmcou(1)=nmcou(1)+1 imcout(nmcou(1),1)=j endif if(ids(j).eq.2)then c c status is missing ordered discrete. c nmcou(2)=nmcou(2)+1 imcout(nmcou(2),2)=j endif if(ids(j).eq.3)then c c status is missing continuous. c nmcou(3)=nmcou(3)+1 imcout(nmcou(3),3)=j endif endif if(imx(i,j).eq.1)then if(ids(j).eq.1)then c c status is observed binary. c nmcou(4)=nmcou(4)+1 imcout(nmcou(4),4)=j endif if(ids(j).eq.2)then c c status is observed ordered discrete. c nmcou(5)=nmcou(5)+1 imcout(nmcou(5),5)=j endif if(ids(j).eq.3)then c c status is observed continuous. c nmcou(6)=nmcou(6)+1 imcout(nmcou(6),6)=j endif endif c c status is always missing. c if(ialmis(j).eq.1)then nmcou(7)=nmcou(7)+1 imcout(nmcou(7),7)=j endif 4 continue c nunob=nmcou(1)+nmcou(2)+nmcou(3) nob=nx-(nunob+nmcou(7)) c if(iopout.eq.2)then write(18,215)nunob 215 format(1x,i3) write(18,214)(imcou(j),j=1,nunob) 214 format(122(1x,i3)) endif c c--------------------------------------------------------------------- c if(nunob.eq.0)then c c no missing variables. update sample moments. c do 90 j=4,6 if(nmcou(j).gt.0)then nmcoutj=nmcou(j) do 88 jj=1,nmcoutj imjt1=imcout(jj,j) if(istan .ne. 1) xtr(imjt1)=x(i,imjt1) if(istan .eq. 1) then if(ids(imjt1).ne.3)xtr(imjt1)=x(i,imjt1) if(ids(imjt1).eq.3)then xu=x(i,imjt1) call undo(xu,imjt1,xut) xtr(imjt1)=xut endif endif 88 continue endif 90 continue do 44 j=4,6 if(nmcou(j).gt.0)then nmcoutj=nmcou(j) do 84 jj=1,nmcoutj imjt1=imcout(jj,j) sob(imjt1,1)=sob(imjt1,1)+1. smean(imjt1,1)=smean(imjt1,1)+xtr(imjt1) if(iouttype.eq.1)then sob(imjt1,2)=sob(imjt1,2)+1. smean(imjt1,2)=smean(imjt1,2)+xtr(imjt1) endif do 45 k=4,6 if(nmcou(k).gt.0)then nmcouk=nmcou(k) do 85 kk=1,nmcouk imkt1=imcout(kk,k) sob2(imjt1,imkt1)= & sob2(imjt1,imkt1)+1. scovt=xtr(imjt1)*xtr(imkt1) if(iouttype.eq.1)then do 94 l=1,3 94 scov(imjt1,imkt1,l)= & scov(imjt1,imkt1,l)+scovt endif if(iouttype.eq.2) & scov(imjt1,imkt1,1)= & scov(imjt1,imkt1,1)+scovt 85 continue endif 45 continue 84 continue endif 44 continue c cx goto 2 goto 97 c endif c c--------------------------------------------------------------------- c c concatenate. c kt=0 do 5 j=1,6 if(nmcou(j).gt.0)then nmcout=nmcou(j) do 6 k=1,nmcout kt=kt+1 6 imcou(kt)=imcout(k,j) endif 5 continue nxadj=kt c c reorder mean and covariance matrix. c do 7 j=1,nxadj ameant(j)=amean(imcou(j)) do 8 k=1,nxadj 8 acovt(j,k)=acov(imcou(j),imcou(k)) 7 continue c c compute mean and covariance matrix of unobserved and observed c discrete latent variables conditional on observed continuous c variables. c npart1=0 do 10 j=1,5 10 npart1=npart1+nmcou(j) npart2=nx-npart1 c if(nmcou(6).eq.0)then c c there are no observed continuous variables. c do 11 j=1,npart1 conmean(j)=ameant(j) do 12 k=1,npart1 12 concov(j,k)=acovt(j,k) 11 continue c endif c if(nmcou(6).gt.0)then c c there are observed continuous variables. invert acovt. c call linv3f(acovt,work1,1,nx,nxs,d1,d2,work2,ier) if(ier.gt.0)then write(6,200)ier,i 200 format(1x,'ier = ',i3,' i = ',i5) stop endif c c partition the inverse of acovt. c do 9 j=1,npart1 conmean(j)=ameant(j) do 13 k=1,npart1 13 concov(j,k)=acovt(j,k) kt=npart1 do 14 k=1,npart2 kt=kt+1 14 cross(j,k)=acovt(j,kt) 9 continue jt=npart1 do 15 j=1,npart2 jt=jt+1 15 dev(j)=x(i,imcou(jt))-ameant(jt) c c invert partitioned covariance matrix to get conditional c covariance matrix. c call linv3f(concov,work1,1,npart1,nxs,d1,d2,work2,ier) if(ier.gt.0)then write(6,200)ier,i stop endif c c adjust conditional mean. c do 16 j=1,npart1 do 17 k=1,npart1 do 18 l=1,npart2 18 conmean(j)=conmean(j)-(concov(j,k)* & cross(k,l)*dev(l)) 17 continue 16 continue c endif c c--------------------------------------------------------------------- c c simulate vector of unobserved variables conditional on observed c variables. c nsuc=0 if(iopout.eq.3)then do 109 j=1,nx nxbnryt(j)=-9 nxorddt(j)=-9 xrealt(j)=-9. if(imx(i,j).eq.1)then if(ids(j).eq.1)nxbnryt(j)=nxbnry(i,j) if(ids(j).eq.2)nxorddt(j)=nxordd(i,j) if(ids(j).eq.3)xrealt(j)=xreal(i,j) endif 109 continue endif do 19 j=1,nunob umean(j)=0. do 20 k=1,nunob 20 ucov(j,k)=0. 19 continue c c--------------------------------------------------------------------- c iseed=1 95 continue c c--------------------------------------------------------------------- c if(iacrj.eq.1)then c c use acceptance-rejection methods to simulate e. c if(nsuc.eq.0)then jk=0 do 22 j=1,npart1 do 23 k=1,j jk=jk+1 23 amat(jk)=concov(j,k) 22 continue call choles(amat,nxs2,nx,cmat) endif 25 continue do 21 j=1,npart1 call nrandom(iseed,ut) iseed=0 21 u(j)=ut call errget(e,npart1,cmat,nxs2,u) do 31 j=1,npart1 31 e(j)=e(j)+conmean(j) c if((nmcou(4)+nmcou(5)).gt.0)then c c check if e satisfies constraints implied by observed discrete c variables. c jt1=nunob if(nmcou(4).gt.0)then nmcout=nmcou(4) do 24 j=1,nmcout jt1=jt1+1 if((x(i,imcou(jt1)).eq.0.).and.(e(j).gt.0.)) & goto 25 if((x(i,imcou(jt1)).eq.1.).and.(e(j).lt.0.)) & goto 25 24 continue endif if(nmcou(5).gt.0)then nmcout=nmcou(5) jt2=nmcou(4) do 26 j=1,nmcout jt1=jt1+1 jt2=jt2+1 imjt1=imcou(jt1) ixt=x(i,imjt1)+0.0001 if(((e(jt2).lt.chat(imjt1,ixt)).and. & (ixt.gt.1)).or.((e(jt2).gt. & chat(imjt1,(ixt+1))).and.((ixt+1).lt. & nchat(imjt1))))goto 25 26 continue endif endif c endif c if(iacrj.eq.2)then ntemp=npart1 jco(1)=0 do 82 j=2,5 82 jco(j)=jco(j-1)+nmcou(j-1) jjj=0 do 65 j=1,5 jt=jtrans(j) nmcoutj=nmcou(jt) if(nmcoutj.gt.0)then jcot=jco(jt) do 83 jj=1,nmcoutj jcot=jcot+1 jjj=jjj+1 amntemp(jjj)=conmean(jcot) kkk=0 do 81 k=1,5 kt=jtrans(k) nmcoutk=nmcou(kt) if(nmcoutk.gt.0)then kcot=jco(kt) do 76 kk=1,nmcoutk kcot=kcot+1 kkk=kkk+1 76 covtemp(jjj,kkk)=concov(jcot,kcot) endif 81 continue 83 continue endif 65 continue c if((nmcou(4)+nmcou(5)).gt.0)then c c use the geweke-hajivassilliou-keane method to simulate e. c ieigchk=1 jt1=nunob if(nmcou(4).gt.0)then nmcout=nmcou(4) do 64 j=1,nmcout if(ieigchk.eq.1)then kl=0 do 57 k=1,ntemp do 57 l=1,k kl=kl+1 57 covsym(kl)=covtemp(k,l) call eigrs(covsym,ntemp,0,eigval,eigvec, & nxs,work1,ier) if(eigval(1).le.0.)then write(6,203)i,ntemp,(eigval(k), & k=1,ntemp) 203 format(1x,'negative eigenvalue: i: ', & i5,3x,'jjj: ',i3,/,122(1x,g15.8)) stop endif endif jt1=jt1+1 ut=rand2(iseed) iseed=0 if(covtemp(1,1).gt.0.)then arg=-amntemp(1)/(covtemp(1,1)**.5) farg=anormc(arg) endif if(covtemp(1,1).le.0.)then if(amntemp(1).lt.0.)farg=0.9999999 if(amntemp(1).gt.0.)farg=0.0000001 if(amntemp(1).eq.0.)farg=0.5 endif ut=ut*farg if(ut.gt..9999999)ut=.9999999 if(ut.lt..0000001)ut=.0000001 call mdnris(ut,arg,ier) if(abs(arg).gt.19.999)then if(arg.gt.19.999)arg=19.999 if(arg.lt.-19.999)arg=-19.999 endif if(x(i,imcou(jt1)).eq.1.)then ut=(ut*(1.-farg))+farg if(ut.gt..9999999)ut=.9999999 if(ut.lt..0000001)ut=.0000001 call mdnris(ut,arg,ier) if(abs(arg).gt.19.999)then if(arg.gt.19.999)arg=19.999 if(arg.lt.-19.999)arg=-19.999 endif endif e(j)=arg call linv3f(covtemp,work1,1,ntemp,nxs,d1,d2, & work2,ier) if(ier.gt.0)then write(6,200)ier,i stop endif do 66 k=2,ntemp 66 crotemp(k)=covtemp(1,k) ntemp=ntemp-1 do 67 k=1,ntemp do 67 l=1,ntemp 67 covtemp(k,l)=covtemp((k+1),(l+1)) call linv3f(covtemp,work1,1,ntemp,nxs,d1,d2, & work2,ier) if(ier.gt.0)then write(6,200)ier,i stop endif devtemp=e(j)-amntemp(1) do 68 k=1,ntemp amntemp(k)=amntemp(k+1) do 98 l=1,ntemp amntemp(k)=amntemp(k)-(covtemp(k,l)* & crotemp(l)*devtemp) 98 continue 68 continue 64 continue endif c if(nmcou(5).gt.0)then nmcout=nmcou(5) do 69 j=1,nmcout jt1=jt1+1 imjt1=imcou(jt1) ixt=x(i,imjt1)+0.0001 ut=rand2(iseed) iseed=0 if(covtemp(1,1).gt.0.)then arg=chat(imjt1,ixt)-amntemp(1)/ & (covtemp(1,1)**.5) farg=anormc(arg) endif if(covtemp(1,1).le.0.)then arg=chat(imjt1,ixt)-amntemp(1) if(arg.lt.0.)farg=0.0000001 if(arg.gt.0.)farg=0.9999999 if(arg.eq.0.)farg=0.5 endif farg1=farg if(covtemp(1,1).gt.0.)then arg=chat(imjt1,(ixt+1))-amntemp(1)/ & (covtemp(1,1)**.5) farg=anormc(arg) endif if(covtemp(1,1).le.0.)then arg=chat(imjt1,(ixt+1))-amntemp(1) if(arg.lt.0.)farg=0.0000001 if(arg.gt.0.)farg=0.9999999 if(arg.eq.0.)farg=0.5 endif farg2=farg ut=(ut*(farg2-farg1))+farg1 if(ut.gt..9999999)ut=.9999999 if(ut.lt..0000001)ut=.0000001 call mdnris(ut,arg,ier) if(abs(arg).gt.19.999)then if(arg.gt.19.999)arg=19.999 if(arg.lt.-19.999)arg=-19.999 endif e(j)=arg call linv3f(covtemp,work1,1,ntemp,nxs,d1,d2, & work2,ier) if(ier.gt.0)then write(6,200)ier,i stop endif do 70 k=2,ntemp 70 crotemp(k)=covtemp(1,k) ntemp=ntemp-1 do 71 k=1,ntemp do 71 l=1,ntemp 71 covtemp(k,l)=covtemp((k+1),(l+1)) call linv3f(covtemp,work1,1,ntemp,nxs,d1,d2, & work2,ier) if(ier.gt.0)then write(6,200)ier,i stop endif devtemp=e(j)-amntemp(1) do 72 k=1,ntemp amntemp(k)=amntemp(k+1) do 72 l=1,ntemp 72 amntemp(k)=amntemp(k)-(covtemp(k,l)* & crotemp(l)*devtemp) 69 continue endif endif do 73 j=1,ntemp call nrandom(iseed,ut) iseed=0 73 u(j)=ut jk=0 do 54 j=1,ntemp do 56 k=1,j jk=jk+1 56 amat(jk)=covtemp(j,k) 54 continue call choles(amat,nxs2,nx,cmat) call errget(et,ntemp,cmat,nxs2,u) do 53 j=1,ntemp 53 e(j)=et(j)+amntemp(j) endif c if(iacrj.eq.0)then jco(1)=0 do 114 j=2,5 114 jco(j)=jco(j-1)+nmcou(j-1) jjj=0 do 115 j=1,5 jt=jtrans(j) nmcoutj=nmcou(jt) if(nmcoutj.gt.0)then jcot=jco(jt) do 116 jj=1,nmcoutj jcot=jcot+1 jjj=jjj+1 amntemp(jjj)=conmean(jcot) kkk=0 do 117 k=1,5 kt=jtrans(k) nmcoutk=nmcou(kt) if(nmcoutk.gt.0)then kcot=jco(kt) do 118 kk=1,nmcoutk kcot=kcot+1 kkk=kkk+1 118 covtemp(jjj,kkk)=concov(jcot,kcot) endif 117 continue 116 continue endif 115 continue c if((nmcou(4)+nmcou(5)).gt.0)then c c use the geweke-hajivassilliou-keane method to simulate e. c ieigchk=1 jt1=nunob if(nmcou(4).gt.0)then nmcout=nmcou(4) do 119 j=1,nmcout if(ieigchk.eq.1)then kl=0 do 120 k=1,j do 120 l=1,k kl=kl+1 120 covsym(kl)=covtemp(k,l) call eigrs(covsym,j,0,eigval,eigvec,nxs, & work1,ier) if(eigval(1).le.0.)then write(6,221)i,j,(eigval(k),k=1,j) 221 format(1x,'negative eigenvalue: i: ', & i5,3x,'j: ',i3,/,122(1x,g15.8)) stop endif endif jt1=jt1+1 ut=rand2(iseed) iseed=0 amnsca=amntemp(j) asigsca=covtemp(j,j) if(j.gt.1)then j1=j-1 do 121 k=1,j1 do 121 l=1,j1 121 covinv(k,l)=covtemp(k,l) d1=-1. call linv3f(covinv,work1,1,j1,nxs,d1,d2, & work2,ier) do 122 k=1,j1 do 122 l=1,j1 amnsca=amnsca+(covtemp(j,k)* & covinv(k,l)*devt(l)) 122 asigsca=asigsca-(covtemp(j,k)* & covinv(k,l)*covtemp(j,l)) endif if(asigsca.gt.0.)then arg=-amnsca/(asigsca**.5) farg=anormc(arg) endif if(asigsca.le.0.)then if(amnsca.lt.0.)farg=0.9999999 if(amnsca.gt.0.)farg=0.0000001 if(amnsca.eq.0.)farg=0.5 endif if(x(i,imcou(jt1)).eq.0.)then ut=ut*farg if(ut.gt..9999999)ut=.9999999 if(ut.lt..0000001)ut=.0000001 call mdnris(ut,arg,ier) if(abs(arg).gt.19.999)then if(arg.gt.19.999)arg=19.999 if(arg.lt.-19.999)arg=-19.999 endif endif if(x(i,imcou(jt1)).eq.1.)then ut=(ut*(1.-farg))+farg if(ut.gt..9999999)ut=.9999999 if(ut.lt..0000001)ut=.0000001 call mdnris(ut,arg,ier) if(abs(arg).gt.19.999)then if(arg.gt.19.999)arg=19.999 if(arg.lt.-19.999)arg=-19.999 endif endif e(j)=((asigsca**.5)*arg)+amnsca devt(j)=e(j)-amntemp(j) 119 continue endif c if(nmcou(5).gt.0)then nmcout=nmcou(5) do 125 j=1,nmcout jt1=jt1+1 imjt1=imcou(jt1) ixt=x(i,imjt1)+0.0001 ut=rand2(iseed) iseed=0 j1=nmcou(4)+j-1 j11=j1+1 amnsca=amntemp(j11) asigsca=covtemp(j11,j11) if((jt1-nunob).gt.1)then do 123 k=1,j1 do 123 l=1,j1 123 covinv(k,l)=covtemp(k,l) d1=-1. call linv3f(covinv,work1,1,j1,nxs,d1,d2, & work2,ier) do 124 k=1,j1 do 124 l=1,j1 amnsca=amnsca+(covtemp(j11,k)* & covinv(k,l)*devt(l)) 124 asigsca=asigsca-(covtemp(j11,k)* & covinv(k,l)*covtemp(j11,l)) endif if(asigsca.gt.0.)then arg=chat(imjt1,ixt)-amnsca/(asigsca**.5) farg=anormc(arg) endif if(asigsca.le.0.)then arg=chat(imjt1,ixt)-amnsca if(arg.lt.0.)farg=0.0000001 if(arg.gt.0.)farg=0.9999999 if(arg.eq.0.)farg=0.5 endif farg1=farg if(asigsca.gt.0.)then arg=chat(imjt1,(ixt+1))-amnsca/ & (asigsca**.5) farg=anormc(arg) endif if(asigsca.le.0.)then arg=chat(imjt1,(ixt+1))-amnsca if(arg.lt.0.)farg=0.0000001 if(arg.gt.0.)farg=0.9999999 if(arg.eq.0.)farg=0.5 endif farg2=farg ut=(ut*(farg2-farg1))+farg1 if(ut.gt..9999999)ut=.9999999 if(ut.lt..0000001)ut=.0000001 call mdnris(ut,arg,ier) if(abs(arg).gt.19.999)then if(arg.gt.19.999)arg=19.999 if(arg.lt.-19.999)arg=-19.999 endif e(j11)=((asigsca**.5)*arg)+amnsca devt(j11)=e(j11)-amntemp(j11) 125 continue endif endif ntemp=nmcou(1)+nmcou(2)+nmcou(3) do 126 j=1,ntemp call nrandom(iseed,ut) iseed=0 126 u(j)=ut jt1=nmcou(4)+nmcou(5) do 127 j=1,ntemp jt1j=jt1+j amnsim(j)=amntemp(jt1j) do 127 k=1,ntemp jt1k=jt1+k 127 asigsim(j,k)=covtemp(jt1j,jt1k) if(jt1.gt.0)then do 128 j=1,jt1 do 128 k=1,jt1 128 covinv(j,k)=covtemp(j,k) d1=-1. call linv3f(covinv,work1,1,jt1,nxs,d1,d2,work2,ier) do 129 j=1,ntemp jt1j=jt1+j do 130 k=1,jt1 do 130 l=1,jt1 130 amnsim(j)=amnsim(j)+(covtemp(jt1j,k)* & covinv(k,l)*devt(l)) do 131 k=1,ntemp jt1k=jt1+k do 132 l=1,jt1 do 132 m=1,jt1 132 asigsim(j,k)=asigsim(j,k)- & (covtemp(jt1j,l)*covinv(l,m)* & covtemp(jt1k,m)) 131 continue 129 continue endif jk=0 do 133 j=1,ntemp do 134 k=1,j jk=jk+1 covsym(jk)=asigsim(j,k) 134 amat(jk)=asigsim(j,k) 133 continue call eigrs(covsym,ntemp,0,eigval,eigvec,nxs,work1,ier) if(eigval(1).le.0.)then write(6,203)i,ntemp,(eigval(k),k=1,ntemp) stop endif call choles(amat,nxs2,ntemp,cmat) call errget(et,ntemp,cmat,nxs2,u) do 135 j=1,ntemp 135 e(j)=et(j)+amnsim(j) endif c c e satisfies the constraints. discretize the discrete unobserved c variables. c if(nmcou(1).gt.0)then nmcout=nmcou(1) do 27 j=1,nmcout if(e(j).le.0.)xt(j)=0. if(e(j).gt.0.)xt(j)=1. 27 continue endif jt=nmcou(1) if(nmcou(2).gt.0)then nmcout=nmcou(2) do 28 j=1,nmcout jt=jt+1 imjt1=imcou(jt) nchatj=nchat(imjt1)-1 do 29 k=1,nchatj if(((e(jt).gt.chat(imjt1,k)).or.(k.eq.1)).and. & ((e(jt).le.chat(imjt1,(k+1))).or.((k+1).eq. & nchat(imjt1))))then xt(jt)=k goto 28 endif 29 continue 28 continue endif c c complete xt. c if(nmcou(3).gt.0)then nmcout=nmcou(3) do 30 j=1,nmcout jt=jt+1 30 xt(jt)=e(jt) endif c c update. c do 89 j=1,nunob imjt1=imcou(j) if(istan .ne. 1) xtr(j)=xt(j) if(istan .eq. 1) then if(ids(imjt1).ne.3)xtr(j)=xt(j) if(ids(imjt1).eq.3)then xu=xt(j) call undo(xu,imjt1,xut) xtr(j)=xut endif endif 89 continue do 32 j=1,nunob umean(j)=umean(j)+xtr(j) do 33 k=1,nunob 33 ucov(j,k)=ucov(j,k)+(xtr(j)*xtr(k)) 32 continue c if(iopout.eq.2)write(18,213)(xtr(j),j=1,nunob) 213 format(122(1x,f6.3)) c c--------------------------------------------------------------------- c 97 continue c c--------------------------------------------------------------------- c c write the new data to an output file. c if ((iopout.eq.3).and. (imonte.eq.1)) then if (nunob .eq. 0) then do 99 j = 1,ndraw 99 write(18,219)i,(x(i,k),k=1,nx) goto 2 endif if (nunob .gt. 0) then do 110 j = 1,nunob imjt1 = imcou(j) x(i,imjt1) = xt(j) if ((ids(imjt1).eq.3).and.(istan.eq.1)) then xu = xt(j) callundo(xu,j,xut) x(i,imjt1) = xut endif 110 continue write(18,219)i,(x(i,k),k=1,nx) endif endif 219 format(1x,i4,15(1x,f8.4)) c c--------------------------------------------------------------------- c nsuc=nsuc+1 if(nsuc.lt.ndraw)goto 95 c c--------------------------------------------------------------------- c c adjust moments. c do 34 j=1,nunob 34 umean(j)=umean(j)/ndraw do 35 j=1,nunob do 36 k=1,nunob 36 ucov(j,k)=(ucov(j,k)/ndraw)-(umean(j)*umean(k)) 35 continue c c update sample moments. c do 37 j=1,nunob imjt1=imcou(j) sob(imjt1,2)=sob(imjt1,2)+1. smean(imjt1,2)=smean(imjt1,2)+umean(j) do 40 k=1,nunob imkt1=imcou(k) scov(imjt1,imkt1,2)=scov(imjt1,imkt1,2)+(umean(j)* & umean(k)) 40 scov(imjt1,imkt1,3)=scov(imjt1,imkt1,3)+ucov(j,k) 37 continue do 86 j=1,nx if(istan .ne. 1) x(i,j)=x(i,j) if(istan .eq. 1) then if(ids(j).eq.3)then xu=x(i,j) call undo(xu,j,xut) x(i,j)=xut endif endif 86 continue c if(nob.gt.0)then jt=nunob do 41 j=1,nob jt=jt+1 imjt1=imcou(jt) sob(imjt1,1)=sob(imjt1,1)+1. smean(imjt1,1)=smean(imjt1,1)+x(i,imjt1) if(iouttype.eq.1)then sob(imjt1,2)=sob(imjt1,2)+1. smean(imjt1,2)=smean(imjt1,2)+x(i,imjt1) do 42 k=1,nunob imkt1=imcou(k) scovt=umean(k)*x(i,imjt1) scov(imjt1,imkt1,2)=scov(imjt1,imkt1,2)+scovt 42 scov(imkt1,imjt1,2)=scov(imkt1,imjt1,2)+scovt endif kt=nunob do 43 k=1,nob kt=kt+1 imkt1=imcou(kt) sob2(imjt1,imkt1)=sob2(imjt1,imkt1)+1. scovt=x(i,imjt1)*x(i,imkt1) scov(imjt1,imkt1,1)=scov(imjt1,imkt1,1)+scovt if(iouttype.eq.1)scov(imjt1,imkt1,2)= & scov(imjt1,imkt1,2)+scovt 43 continue 41 continue c endif c endif c if(iopout.eq.1)then write(6,211)i 211 format(1x,'imputed values for observation ',i4,/1x, & 'var no',2x,'type',6x,'mean',6x,'std. dev.') do 87 j=1,nunob imjt1=imcou(j) if(ucov(j,j).gt.0.)sucov=ucov(j,j)**.5 if(ucov(j,j).le.0.)sucov=0. 87 write(6,212)imjt1,aids(ids(imjt1)),umean(j),sucov 212 format(1x,i3,5x,a8,2(2x,f8.5)) c endif c if(iopout.eq.0)then if(((i/10)*10).eq.i)write(6,216)i 216 format(1x,'i = ',i5) endif c 2 continue c c adjust sample moments. c do 46 i=1,nx if(ialmis(i).eq.0)then smean(i,1)=smean(i,1)/sob(i,1) smean(i,2)=smean(i,2)/sob(i,2) endif 46 continue do 47 i=1,nx if(ialmis(i).eq.0)then do 48 j=1,nx if(ialmis(j).eq.0)then scov(i,j,1)=(scov(i,j,1)/sob2(i,j))-(smean(i,1)* & smean(j,1)) if(iouttype.eq.1)scov(i,j,2)=(scov(i,j,2)/ & sob(i,2))-(smean(i,2)*smean(j,2)) if((iouttype.eq.2).and.(i.eq.j))then if(sob(i,2).gt.0.)scov(i,i,2)=(scov(i,i,2)/ & sob(i,2))-(smean(i,2)*smean(i,2)) if(sob(i,2).le.0.)scov(i,i,2)=0. endif if(iouttype.eq.1)then if(sob(i,2).gt.sob(i,1))scov(i,j,3)=scov(i,j,3)/ & (sob(i,2)-sob(i,1)) if(sob(i,2).le.sob(i,1))scov(i,j,3)=0. endif if((iouttype.eq.2).and.(i.eq.j))then if(sob(i,2).gt.0.)scov(i,i,3)=scov(i,i,3)/ & sob(i,2) if(sob(i,2).le.0.)scov(i,i,3)=0. endif endif 48 continue endif 47 continue c c output. c if(iouttype.eq.1)write(15,218)ii 218 format(1x,'means and standard deviations for type ',i2,/,9x, & 'observed sample',6x,'full sample',10x,'std dev of mean',/,1x, & 'var no',2x,'mean',6x,'std dev') if(iouttype.eq.2)write(15,201)ii 201 format(1x,'means and standard deviations for type ',i2,/,9x, & 'observed sample',6x,'missing sample',7x,'std dev of mean',/, & 1x,'var no',2x,'mean',6x,'std dev') do 49 i=1,nx if(ialmis(i).eq.0)then do 55 j=1,3 55 std(j)=scov(i,i,j)**.5 write(15,217)i,(smean(i,j),std(j),j=1,2),std(3), & (sob(i,j),j=1,2) endif 49 continue 217 format(1x,i3,5x,f8.4,2x,f8.4,4x,f8.4,2x,f8.4,4x,f8.4,/,9x, & f8.1,14x,f8.1) 52 continue return end c c----------------------------------------------------------------------- c c imsl routine name - linv3f c c----------------------------------------------------------------------- c c computer - prim77/single c c latest revision - june 1, 1982 c c purpose - in place inverse, equation solution, and/or c determinant evaluation - full storage mode c c usage - call linv3f (a,b,ijob,n,ia,d1,d2,wkarea,ier) c c arguments a - input/output matrix of dimension n by n. see c parameter ijob. c b - input/output vector of length n when ijob = c 2 or 3. otherwise, b is not used. c on input, b contains the right hand side of c of the equation ax = b. c on output, the solution x replaces b. c ijob - input option parameter. ijob = i implies: c i = 1, invert matrix a. a is replaced by c its inverse. c i = 2, solve the equation ax = b. a is c replaced by the lu decomposition of a c rowwise permutation of a, where u is c upper triangular and l is lower c triangular with unit diagonal. c the unit diagonal of l is not stored. c i = 3, solve ax = b and invert matrix a. c a is replaced by its inverse. c i = 4, compute the determinant of a. c a is replaced by the lu decomposition c of a rowwise permutation of a. c n - order of a. (input) c ia - row dimension of matrix a exactly as c specified in the dimension statement in the c calling program. (input) c d1 - input/output. if the d1 and d2 components c d2 of determinant(a) = d1*2**d2 are desired, c input d1.ge.0. otherwise, input d1.lt.0. c d2 is never input. c wkarea - work area of length at least 2*n for ijob=1 c or ijob=3. c work area of length at least n for ijob=2 c or ijob=4. c ier - error parameter. (output) c warning with fix c ier = 65 indicates that ijob was less than c 1 or greater than 4. ijob is assumed to c be 4. c terminal error c ier = 130 indicates that matrix a is c algorithmically singular. (see the c chapter l prelude). c c precision/hardware - single and double/h32 c - single/h36,h48,h60 c c reqd. imsl routines - ludatn,luelmn,uertst,ugetio c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1978 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c----------------------------------------------------------------------- c subroutine linv3f (a,b,ijob,n,ia,d1,d2,wkarea,ier) c implicit real*8 (a-h,o-z) c implicit real (a-h,o-z) real*8 a(ia,n),b(n),wkarea(300) c real*4 a(ia,n),b(n),wkarea(300) data zero/0.d0/,one/1.d0/ c first executable statement c lu decomposition of a call ludatn (a,ia,n,a,ia,0,c1,c2,wkarea,wkarea,wa,ier) if (d1 .lt. zero .and. ijob .ge. 1 .and. ijob .lt. 4) go to 5 d1 = c1 d2 = c2 5 if (ier .ge. 128) go to 60 if (ijob .le. 0 .or. ijob .gt. 4) go to 55 c solve ax = b if (ijob .eq. 2 .or . ijob .eq. 3) call luelmn (a,ia,n,b,wkarea,b) if (ijob .ne. 1 .and. ijob .ne. 3) go to 9005 c matrix inversion a(n,n) = one/a(n,n) nm1 = n-1 if (nm1 .lt. 1) go to 9005 do 40 ii=1,nm1 l = n-ii m = l+1 do 15 i=m,n sum = zero do 10 k=m,n sum = sum-a(i,k)*a(k,l) 10 continue wkarea(n+i) = sum 15 continue do 20 i=m,n a(i,l) = wkarea(n+i) 20 continue do 30 j=l,n sum = zero if (j .eq. l) sum = one do 25 k=m,n sum = sum-a(l,k)*a(k,j) 25 continue wkarea(n+j) = sum/a(l,l) 30 continue do 35 j=l,n a(l,j) = wkarea(n+j) 35 continue 40 continue c permute columns of a inverse do 50 i=1,n j = n-i+1 jp = wkarea(j) if (j .eq. jp) go to 50 do 45 k=1,n c = a(k,jp) a(k,jp) = a(k,j) a(k,j) = c 45 continue 50 continue go to 9005 55 continue c warning with fix - ijob was set c incorrectly ier = 65 go to 9000 c terminal error - matrix a is c algorithmically singular 60 ier = 130 9000 continue c call uertst(ier,'linv3f') 9005 return end c c--------------------------------------------------------------------- c c imsl routine name - ludatn c c----------------------------------------------------------------------- c c computer - prim77/single c c latest revision - june 1, 1982 c c purpose - nucleus called only by imsl subroutine leqt2f c c precision/hardware - single and double/h32 c - single/h36,h48,h60 c c reqd. imsl routines - uertst,ugetio c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1982 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c----------------------------------------------------------------------- c subroutine ludatn (a,ia,n,xlu,ilu,idgt,d1,d2,apvt,equil,wa,ier) implicit real*8 (a-h,o-z) c implicit real*4 (a-h,o-z) c real*8 a(ia,n),xlu(ilu,ilu),apvt(ilu),equil(n) c real*4 a(ia,n),xlu(ilu,ilu),apvt(ilu),equil(n) data zero,one,four,sixtn,sixth/0.d0,1.d0,4.d0,16.d0,.0625d0/ c first executable statement c initialization ier = 0 rn = n wrel = zero d1 = one d2 = zero biga = zero do 10 i=1,n big = zero do 5 j=1,n p = a(i,j) xlu(i,j) = p c p = dabs(p) p = abs(p) if (p .gt. big) big = p 5 continue if (big .gt. biga) biga = big if (big .eq. zero) go to 110 equil(i) = one/big 10 continue do 105 j=1,n jm1 = j-1 if (jm1 .lt. 1) go to 40 c compute u(i,j), i=1,...,j-1 do 35 i=1,jm1 sum = xlu(i,j) im1 = i-1 if (idgt .eq. 0) go to 25 c with accuracy test c ai = dabs(sum) ai = abs(sum) wi = zero if (im1 .lt. 1) go to 20 do 15 k=1,im1 t = xlu(i,k)*xlu(k,j) sum = sum-t c wi = wi+dabs(t) wi = wi+abs(t) 15 continue xlu(i,j) = sum c 20 wi = wi+dabs(sum) 20 wi = wi+abs(sum) if (ai .eq. zero) ai = biga test = wi/ai if (test .gt. wrel) wrel = test go to 35 c without accuracy 25 if (im1 .lt. 1) go to 35 do 30 k=1,im1 sum = sum-xlu(i,k)*xlu(k,j) 30 continue xlu(i,j) = sum 35 continue 40 p = zero c compute u(j,j) and l(i,j), i=j+1,..., do 70 i=j,n sum = xlu(i,j) if (idgt .eq. 0) go to 55 c with accuracy test c ai = dabs(sum) ai = abs(sum) wi = zero if (jm1 .lt. 1) go to 50 do 45 k=1,jm1 t = xlu(i,k)*xlu(k,j) sum = sum-t c wi = wi+dabs(t) wi = wi+abs(t) 45 continue xlu(i,j) = sum c 50 wi = wi+dabs(sum) 50 wi = wi+abs(sum) if (ai .eq. zero) ai = biga test = wi/ai if (test .gt. wrel) wrel = test go to 65 c without accuracy test 55 if (jm1 .lt. 1) go to 65 do 60 k=1,jm1 sum = sum-xlu(i,k)*xlu(k,j) 60 continue xlu(i,j) = sum c 65 q = equil(i)*dabs(sum) 65 q = equil(i)*abs(sum) if (p .ge. q) go to 70 p = q imax = i 70 continue c test for algorithmic singularity q = rn+p 71 if (q .eq. rn) go to 110 if (j .eq. imax) go to 80 c interchange rows j and imax d1 = -d1 do 75 k=1,n p = xlu(imax,k) xlu(imax,k) = xlu(j,k) xlu(j,k) = p 75 continue equil(imax) = equil(j) 80 apvt(j) = imax d1 = d1*xlu(j,j) c 85 if (dabs(d1) .le. one) go to 90 85 if (abs(d1) .le. one) go to 90 d1 = d1*sixth d2 = d2+four go to 85 c 90 if (dabs(d1) .ge. sixth) go to 95 90 if (abs(d1) .ge. sixth) go to 95 d1 = d1*sixtn d2 = d2-four go to 90 95 continue jp1 = j+1 if (jp1 .gt. n) go to 105 c divide by pivot element u(j,j) p = xlu(j,j) do 100 i=jp1,n xlu(i,j) = xlu(i,j)/p 100 continue 105 continue c perform accuracy test if (idgt .eq. 0) go to 9005 p = 3*n+3 wa = p*wrel q = wa+10.0**(-idgt) 106 if (q .ne. wa) go to 9005 ier = 34 go to 9000 c algorithmic singularity 110 ier = 129 d1 = zero d2 = zero 9000 continue c print error c call uertst(ier,'ludatn') 9005 return end c c--------------------------------------------------------------------- c c imsl routine name - luelmn c c----------------------------------------------------------------------- c c computer - prim77/single c c latest revision - june 1, 1982 c c purpose - nucleus called only by imsl subroutine leqt2f c c reqd. imsl routines - none required c c notation - information on special notation and c conventions is available in the manual c introduction or through imsl routine uhelp c c copyright - 1982 by imsl, inc. all rights reserved. c c warranty - imsl warrants only that imsl testing has been c applied to this code. no other warranty, c expressed or implied, is applicable. c c----------------------------------------------------------------------- c subroutine luelmn (a,ia,n,b,apvt,x) c implicit real*8 (a-h,o-z) implicit real*4 (a-h,o-z) c real*8 a(ia,n),b(ia),apvt(ia),x(ia) c real*4 a(ia,n),b(ia),apvt(ia),x(ia) c first executable statement c solve ly = b for y do 5 i=1,n 5 x(i) = b(i) iw = 0 do 20 i=1,n ip = apvt(i) sum = x(ip) x(ip) = x(i) if (iw .eq. 0) go to 15 im1 = i-1 do 10 j=iw,im1 sum = sum-a(i,j)*x(j) 10 continue go to 20 15 if (sum .ne. 0.) iw = i 20 x(i) = sum c solve ux = y for x do 30 ib=1,n i = n+1-ib ip1 = i+1 sum = x(i) if (ip1 .gt. n) go to 30 do 25 j=ip1,n sum = sum-a(i,j)*x(j) 25 continue 30 x(i) = sum/a(i,i) return end c c--------------------------------------------------------------------- c function rand2(idum) dimension v(97) data iseed/0/ if(idum.ne.0) then temp=rand() do 1 j=1,97 temp=rand() 1 continue do 2 k=1,97 v(k)=rand() 2 continue y=rand() endif j=1+int(97*y) y=v(j) rand2=y v(j)=rand() return end c c--------------------------------------------------------------------- c