c new version stripped, but multi-group c only qu.newton. no toeters en bellen c added april 03 EN(tau) c npsol c c missing data version c eske version aug 06 c program mmdx12 implicit double precision (a-h,o-z),integer (i-n) c parameter (maxg=2,maxn=10000,maxv=12,maxc=1,maxe=12,mpar=25) parameter (ncnln=0,nrowj=1) c parameter (mvt=maxv+maxc+1) parameter (nclin=0,nctotl=(mpar+nclin+ncnln)) parameter (nrowa=1,liwork=3*mpar+nclin+2*ncnln,nhmax=mpar+ncnln) parameter (lwork=2*mpar*mpar+(mpar*nclin)*(2*mpar*ncnln)+ & 20*mpar+11*nclin+21*ncnln) c npsol double precision a(nrowa,mpar),bl(nctotl),bu(nctotl), & cjac(nrowj,mpar),clamda(nctotl),work(lwork),cn(nrowj) integer istate(nctotl),iwork(liwork),inform,iter double precision x(mpar),g(mpar) c hessian double precision hess(nhmax,nhmax) double precision countm(maxg,maxv) c ----------------------------------------------------------------- double precision objf,zero,one,half,shalf,two,pi,crit,oldobjf double precision ctol,ftol,tmp,fml1,objem,entrop,xci(4),se,pi2 double precision ent character version*27,inf1*12,outf*12,inf2*12,oinf2*12 character name(6)*2,tag*2 c ----------------------------------------------------------------- c common / summ / sigma(maxg,maxv,maxv),xmean(maxg,maxv) c common / arra1 / obs(maxn,mvt),xmis,wmis(maxv,maxv) common / arra2 / cov(maxg,maxc,maxv,maxv),xm(maxg,maxc,maxv) common / iarra1 / misindi(maxn),misgr(maxg) c common / wrkm / wmat1(maxv,maxv),wmat2(maxv,maxv), & wmat3(maxv,maxv),wmat4(maxv,maxv) common / wrkv / wvec1(maxv),wvec2(maxv),wvec3(maxv) c common / xlis / xmat(maxg,maxc,6,maxv,maxv),p(maxg,maxc) common / ilis / imat(maxg,maxc,6,maxv,maxv), & idim(maxg,6,2),ip(maxg,maxc),ipat(maxg,6) c common / scala / np,ng,npg(maxg),npc(0:maxg), & nv(maxg),nc(maxg),ne(maxg),itqn common / bits / pi2,half,shalf,one,zero,two common / names / name c ----------------------------------------------------------------- external objfun,fdh,writem,npsol,confun,getpp external readm,fdh2,invert,cholinv,lastwo external hessing,gethes,infofit,teltime,check1,dista c ----------------------------------------------------------------- data zero,half,shalf,one,two,pi & / 0.d0,0.5d0,-0.5d0,1.d0,2.d0,3.141592653589793d0 / data ctol,dcrit / 0.d0,1.d-2 / data name /'LY','BE','PS','TE','NU','AL'/ data version /'VERSION APRIL/02 MG NPSOL '/ c ----------------------------------------------------------------- itqn=0 open(33,file='npsol.in') open(9,file='npsol.out') c ----------------------------------------------------------------- open(17,file='scrax') write(*,993) maxg,maxn,maxv,maxc,maxe,mpar,version c.............................................i/o with user 494 WRITE(*,'(2x,A\)') ' NAME INPUT ' READ(*,'(A12)') INF1 OPEN(1,FILE=INF1,STATUS='OLD',ERR=494) REWIND 1 480 WRITE(*,'(2x,A\)') ' NAME OUTPUT ' READ(*,'(A12)') outf OPEN(16,FILE=outf,status='new',err=491) goto 492 491 write(*,'(2x,a\)') ' ! OVERWRITE OUTFILE (y/n) ? ' read(*,'(a1)') tag if (tag.ne.'y') goto 480 open(16,file=outf) 492 REWIND 16 c .............................................end of conversation call npfile(33,inform) if (inform.ne.0) write(16,*) ' error reading npsol.in ' c open(4,file='postp.dat') c time of day ..... ce call teltime(' start') ce write(16,993) maxg,maxn,maxv,maxc,maxe,mpar,version c write(*,'(2x,a\)') ' number of reps number of sels ' read(*,*) nrep , nsel c - - - - - - - - - - - - - - - do 5000 isel=1,nsel do 5000 irep=1,nrep write(*,*) isel, irep rewind 1 ig=1 np=0 npc(0)=0 do 18 i=1,maxg npc(i)=0 misgr(i)=0 do 18 j=1,maxv countm(i,j)=0.d0 18 continue tmp=0.d0 c read first four lines in group 1 400 continue ce400 write(16,994) ig 994 format('INPUT: GROUP ',i2,/) do 4000 ig=1,2 call reader(ig,inf2,npg(ig),nv(ig),nc(ig),ne(ig), & ng,iwk,maxqn,ipat,xmis,maxg) ftol=one/(10.d0**iwk) if (ig.eq.1.and.irep.eq.1) open(2,file=inf2,status='old',err=516) if (ig.eq.2.and.irep.eq.1) open(3,file=inf2,status='old',err=516) np=np+npg(ig) do 2 j=1,ig npc(ig)=npc(ig)+npg(j) 2 continue do 3 i=npc(ig-1)+1,npc(ig) if (ig.eq.1) read(2,*,err=515)(obs(i,j),j=1,nv(ig)) if (ig.eq.2) read(3,*,err=515)(obs(i,j),j=1,nv(ig)) misindi(i)=0 do 10 j=1,nv(ig) if (obs(i,j).eq.xmis) then countm(ig,j)=countm(ig,j)+1 misindi(i)=misindi(i)+1 misgr(ig)=1 endif 10 continue 3 continue c soms checks --------------------------------------------------- 1 call check1(maxqn,iwk,ne(ig),maxe,nc(ig), & maxn,npg(ig),maxc,nv(ig),maxv,ipat,name,ng,maxg) c write findings in first 4 lines ce write(16,991) ig,npg(ig),nv(ig),ne(ig),nc(ig), ce & (ipat(ig,i),i=1,6),iwk,maxqn ifail=0 c read matrices and proportions ------------------------------------ call readm(ig,nv(ig),ne(ig),nc,idim,ipat,xmat,imat,ng,maxg, & maxc,maxv,p,ip,iwork,liwork,name,npar,mpar,ifail) c some checks ------------------------------------------------------ 2 if (ifail.eq.-9) stop ' unknown matrix name...see output' if (npar.gt.mpar) stop ' npar .gt. mpar (setpar.inc) ' c ce write(16,665) inf1,outf,inf2 c c if (ig.lt.ng) then c ig=ig+1 c goto 400 c endif 4000 continue c c write ce write(16,932) npar if (ncnln.ne.0) write(16,*) 'NUMBER OF NONLIN CONSTRAINTS ',ncnln ce call writem(1) c three scalars pi2=two*pi ntot=npar+ncnln c read start values + bounds if (npar.gt.0) then read(1,*,err=512) (x(i),i=1,npar) read(1,*,err=513) (bl(i),i=1,ntot) read(1,*,err=514) (bu(i),i=1,ntot) endif c read raw data free format c write(16,*) c write(16,*) 'STARTING VALUES AND FIXED VALUES ' c c call objfun(-200,npar,x,objf,g,nstate) c call writem(0) c if (npar.eq.0) then call objfun(-600,npar,x,objf,g,nstate) write(16,967) -objf stop ' no free parameters ! ' endif c--------------------------------------------------------- c ---> check 2 insert here from check.f GRAD CHECK loglikelihood c write(16,*) ' grad check 1:' c call fdh2(x,g,work,npar) c do 603 i=1,npar c tmp=(g(i)-work(i)) c603 write(16,992) i,g(i),work(i),tmp c.......method method..................................0 0 0 0 0 c ce write(16,941) itqn=0 ce call teltime(' start q-n') ce call npsol(npar,nclin,ncnln,nrowa,nrowj,nhmax,a,bl,bu, & confun,objfun,inform,iter,istate,cn,cjac,clamda, & objf,g,hess,x,iwork,liwork,work,lwork) write(16,'(1x,i3,1x,i6)') irep,isel do 19 ig=1,ng 19 write(16,928) ((countm(ig,j)/npg(ig)),j=1,nv(ig)) 928 format(10(f6.3,1x)) write(16,977) inform,2.d0*objf write(16,987) (x(i),i=1,npar) 977 format(1x,i6,1x,f12.3) 5000 continue stop ' OK ? Definitely maybe ! Bye bye... ' c errors 512 write(16,*) ' read error reading starting values ' call lastwo 513 write(16,*) ' read error reading lower bounds ' call lastwo 514 write(16,*) ' read error reading upper bounds ' call lastwo 515 write(16,*) ' error in reading group/case ',ig,i stop ' error in reading data file. see output file ' 516 write(16,*) ' cannot find data file specified as FI=',inf2 stop ' cannot find data file ' c FORMATS 995 format(/,'Q-NEWTON MAXIMIZATION OF LOGL. OF MIXTURE MODEL') 942 format(13(f5.2,1x)) 992 format('P ',i2,' EXACT G ',f15.7,' FD G ',f15.7,' DEV ',f10.6) 932 format(/,'NUMBER OF PARAMETERS ',i2) 881 format( & '-> PAR ',i2,' EST ',f10.5,' GRAD ',f10.5,' > CRIT=',f5.3) 941 format(/,30('-'),' RESULTS ',30('-'),/) 991 format(/,'MODEL: ',/,' - GROUP ',i2,/, &' - NP=',i4,' NV=',i2,' NE=',i2,' NC=',i2,/ &' - PAT LY ',i1,' BE ',i1,' PS ',i1,' TE ',i1, &' NU ',i1,' AL ',i1,/'CONVERGENCE CRIT:',/, &' - QN CRIT 10^-',i2,' QN MAXIT ',i4) 665 format(/, &'INPUT: ',a10,' OUTPUT: ',a10,' DATA: ',a10,/) 948 format('CONVERGED AT IT ',i3,/,' CHI^2 (LISREL) ',f8.2, &' LOGL (EM) ',f10.2) 993 format(1x,74('-'),/, &' CURRENT SETTINGS: NG ',i2,' NP ',i6,' NV ', & i2,' NC ',i1,' NE ',i2,' NPAR ',i3, &/,1x,74('-'),/, &' MULTI-GROUP MULTIVARIATE NORMAL MIXTURE ',/, &' SUBJECT TO STRUCTURAL EQUATION MODELING ',/,1x,a27,/, &1x,74('-'),/) 999 format(7(1x,f10.4)) 998 format(5(1x,i10)) 987 format(3(7(1x,f10.4)/)) 967 format(/,'LOGL FOR FIXED PARAMETERS ',f10.3) c977 format('LOGL ',f10.3,' ML ESTIMATES (Q-N):') 978 format &('LOGL ',f10.3,' CHI^2 (EM BASELINE) ',f10.3,' ML ESTS (Q-N):') 957 format('LOGL ',f10.3,' ML ESTIMATES (EM):') 983 format('GROUP ',i2,' EN(tau) ',f8.4, &' ENTROPY-BASED MEASURE OF FUZZINESS ',f6.3) 913 format(/,'REQUESTED 95% CI ON PARAMETER(S) ',10(i4,1x),/) 929 format &(/,'GROUP ',i2, &' * WARNING: NO ENTROPY-BASED MEASURE..ZERO POSTP.PROBS...',/) end