Program g12zafe ! G12ZAF Example Program Text. ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g11caf, g12zaf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: dev, tol Integer :: cm, i, ifail, ip, iprint, ldz, lisi, & lwk, m, maxit, mxn, n, nd, ns, num, & nxs ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: b(:), cov(:), sc(:), se(:), t(:), & tp(:), wk(:), x(:,:), z(:,:) Integer, Allocatable :: cnt(:), ic(:), id(:), irs(:), & isi(:), isz(:), ixs(:), nca(:), nct(:) ! .. Intrinsic Procedures .. Intrinsic :: count, maxval ! .. Executable Statements .. Write (nout,*) 'G12ZAF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) ! Read in problem size Read (nin,*) n, m, ns, maxit, iprint If (ns>0) Then lisi = n Else lisi = 0 End If ldz = n Allocate (z(ldz,m),isz(m),t(n),ic(n),isi(lisi),tp(n),irs(n)) ! Read in the data If (ns>0) Then Read (nin,*)(t(i),z(i,1:m),ic(i),isi(i),i=1,n) Else Read (nin,*)(t(i),z(i,1:m),ic(i),i=1,n) End If ! Read in the variable indicator Read (nin,*) isz(1:m) ! Calculate number of parameters in the model ip = count(isz(1:m)>0) ! Call the routine once to calculate size of MXN ... ! Dummy allocation mxn = 0 Allocate (x(mxn,ip),id(mxn),ixs(mxn)) ! Call G12ZAF to calculate MXN ifail = 1 Call g12zaf(n,m,ns,z,ldz,isz,ip,t,ic,isi,num,ixs,nxs,x,mxn,id,nd,tp,irs, & ifail) If (ifail/=0 .And. ifail/=3) Then Go To 100 End If ! Required size for MXN is returned in NUM, so reallocate memory mxn = num Deallocate (x,id,ixs) Allocate (x(mxn,ip),id(mxn),ixs(mxn)) ! Create risk set ifail = 0 Call g12zaf(n,m,ns,z,ldz,isz,ip,t,ic,isi,num,ixs,nxs,x,mxn,id,nd,tp,irs, & ifail) Allocate (cnt(nxs),b(ip),se(ip),sc(ip),nca(nxs),nct(nxs),cov(ip*(ip+ & 1)/2)) ! Set tolerance tol = 1.0E-5_nag_wp ! Read in initial parameter estimates Read (nin,*) b(1:ip) ! Count the number of observations in each stratum cnt(1:nxs) = 0 Do i = 1, num cnt(ixs(i)) = cnt(ixs(i)) + 1 End Do cm = maxval(cnt(1:nxs)) lwk = ip*num + (cm+1)*(ip+1)*(ip+2)/2 + cm Allocate (wk(lwk)) ! Get parameter estimates from conditional logistic analysis ifail = 0 Call g11caf(num,ip,nxs,x,mxn,isz,ip,id,ixs,dev,b,se,sc,cov,nca,nct,tol, & maxit,iprint,wk,lwk,ifail) ! Display results Write (nout,*) ' Parameter Estimate', ' Standard Error' Write (nout,*) Write (nout,99999)(i,b(i),se(i),i=1,ip) 100 Continue 99999 Format (I6,10X,F8.4,10X,F8.4) End Program g12zafe