c Simulation study c Generate IOMTM(2) for ordinal data c arbitrary time*beta =0.1d0,-1.5d0,1.2d0,-0.4d0,0.8d0,-0.3d0,-1.d0 c Fit IOMTM(1) c using linear and quadratic c Simulated data $Debug use msimsl parameter (IT=8,K1=4,ialph=2,norder=2,M=1000,M1=5000, & N=500,nx=10,nalpha=1,nalpha1=1) integer nr(N),id(K1,IT,N),r(K1,IT,N),ir(1,K1),idrop(1),t double precision b0(K1-1+nx-1),ba0((K1-1)*nalpha*ialph*norder), & b1(K1-1+nx-1),ba1((K1-1)*nalpha*ialph*norder), & ba11((K1-1)*nalpha1*ialph), & x(nx,IT,N),z(nalpha,IT,N),yv(K1),y(IT,N),beta1(nx,K1), & gamz(ialph,norder,K1-1,IT,N),pc(K1),pm(K1,IT,N),pi(K1,IT,N), & ralpha0(nalpha,ialph,K1-1,norder),triangle(K1-1,IT,N), & pj(K1,K1,IT,N),y1(IT,N),pdrop(2) data b0 / -1.d0,0.7d0,2.d0, & 0.1d0,-1.5d0,1.2d0,-0.4d0,0.8d0,-0.3d0,-1.d0, & 0.1d0,-0.5d0 / data ba11 / -0.1d0, 0.20d0, 0.30d0, & -0.2d0,-0.70d0, -0.30d0 / c data ba11 / -1.8d0,-1.1d0,-0.47d0, c & 0.07d0,-0.4d0,-0.47d0, c & -1.8d0,-1.1d0,-0.47d0, c & 0.07d0,-0.4d0,-0.47d0 / c data ba0 / -0.1d0,-0.3d0,0.5d0,-0.1d0,0.1d0,0.2d0, c & -0.1d0,-0.1d0,0.3d0,-0.2d0,0.2d0,0.1d0, c & -0.1d0,-0.2d0,0.5d0,-0.1d0,0.3d0,0.1d0, c & -0.1d0,-0.4d0,0.3d0,-0.2d0,0.2d0,0.2d0 / c data ba0 / -0.1d0, 0.20d0, 0.3d0, 0.1d0,-0.20d0, 0.5d0, c & -0.2d0,-0.70d0,-0.3d0,-0.1d0, 0.20d0, 0.1d0 / data ba0 / -0.1d0, 0.1d0, 0.2d0, -0.2d0, 0.1d0, 0.1d0, & -0.5d0,-0.70d0,-0.5d0,-0.4d0, 0.50d0, 0.6d0 / data yv /-0.9d0,-0.3d0,0.3d0,0.9d0/ do 20 i=1,N do 23 t=1,IT c constant term x(1,t,i)=1.d0 c Time do 19 ig=1,IT-1 x(1+ig,t,i)=0.d0 if(t.eq.ig) then x(1+ig,t,i)=1.d0 endif 19 continue c write(*,*)(x(ik,t,i),ik=1,IT) c read(*,*)iii c Arm transformation (IFL=control group, FOLFOX, IROX=treatment group) if(i.lt.(N/3))then c group=1 x(IT+1,t,i)=0.d0 x(IT+2,t,i)=0.d0 else if(i.lt.(N/3*2))then c group=2 x(IT+1,t,i)=1.d0 x(IT+2,t,i)=0.d0 else c group=3 x(IT+1,t,i)=0.d0 x(IT+2,t,i)=1.d0 endif 23 continue 20 continue open(3,file='Fit-MTM1_Sim-MTM2_MAR.out') do 70 iter=1,M write(*,*)'iter=',iter do 24 ig=1,K1-1+nx-1 b1(ig)=b0(ig) 24 continue do 26 ig=1,(K1-1)*nalpha*ialph*norder ba1(ig)=ba0(ig) 26 continue do 25 i=1,N nr(i)=IT do 25 t=1,nr(i) z(1,t,i)=1.d0 c z(2,t,i)=x(2,t,i) c z(3,t,i)=x(3,t,i) c z(4,t,i)=x(4,t,i) 25 continue do 27 l=1,norder do 27 j=1,ialph do 27 k=1,K1-1 do 27 ig=1,nalpha ralpha0(ig,j,k,l)= & ba1((K1-1)*ialph*nalpha*(l-1)+ialph*nalpha*(k-1)+ & (j-1)*nalpha+ig) 27 continue do 32 k=1,K1 if(k.lt.K1) then beta1(1,k)=b0(k) do 34 ig=1,nx-1 beta1(ig+1,k)=b1(K1-1+ig) 34 continue else do 36 ig=1,nx beta1(ig,K1)=0.d0 36 continue endif 32 continue c Calculate P^M call calpm(beta1,x,N,nr,IT,K1,nx,pm) c Calculate pi call calpi(pm,N,nr,IT,K1,pi) c Calculate \triangle_{itk} given \beta_0^(n-1), \beta^(n-1) call caltri1(yv,z,pi,pj,ralpha0,gamz,N,nr,IT,K1,ialph,nalpha, & norder,triangle,2) c Calculate P_{itk}^c do 40 i=1,N do ig=1,K1 pc(ig)=pi(ig,1,i) enddo CALL RNMTN (1,1,K1,real(pc),IR,1) do 41 k=1,K1 is=0 do ig1=1,k id(ig1,1,i)=ir(1,ig1) is=is+ir(1,ig1) enddo r(k,1,i)=is if(ir(1,k).eq.1) then y(1,i)=yv(k) y1(1,i)=dble(k)-1.d0 endif 41 continue if(nr(i).ge.2) then do 42 t=2,nr(i) if(t.ge.3) then pdrop(1)=dexp(-2.5d0+0.6d0*y1(t-1,i)+0.3d0*y1(t-2,i))/ & (1.d0+dexp(-2.5d0+0.6d0*y1(t-1,i)+0.3d0*y1(t-2,i))) pdrop(2)=1.d0/ & (1.d0+dexp(-2.5d0+0.6d0*y1(t-1,i)+0.3d0*y1(t-2,i))) call RNBIN(1,1,real(pdrop(1)),idrop) if(idrop(1).eq.1) then nr(i)=t-1 goto 40 endif endif if(t.eq.2) then sumexp=0.d0 do 44 l=1,K1-1 sumexp=sumexp+dexp(triangle(l,t,i)+ & gamz(1,1,l,t,i)*y(t-1,i)+gamz(2,1,l,t,i)*y(t-1,i)**2.d0) 44 continue do 46 k=1,K1-1 pc(k)=dexp(triangle(k,t,i)+ & gamz(1,1,k,t,i)*y(t-1,i)+gamz(2,1,k,t,i)*y(t-1,i)**2.d0)/ & (1.d0+sumexp) 46 continue pc(K1)=1.d0/(1.d0+sumexp) c write(*,*)i,t,pc CALL RNMTN (1,1,K1,real(pc),IR,1) do 47 k=1,K1 is=0 do ig1=1,k id(ig1,t,i)=ir(1,ig1) is=is+ir(1,ig1) enddo r(k,t,i)=is if(ir(1,k).eq.1) then y(t,i)=yv(k) y1(t,i)=dble(k)-1.d0 endif 47 continue else sumexp=0.d0 do 50 l=1,K1-1 sumexp=sumexp+dexp(triangle(l,t,i)+ & gamz(1,1,l,t,i)*y(t-1,i)+gamz(2,1,l,t,i)*y(t-1,i)**2.d0+ & gamz(1,2,l,t,i)*y(t-2,i)+gamz(2,2,l,t,i)*y(t-2,i)**2.d0) 50 continue do 52 k=1,K1-1 pc(k)=dexp(triangle(k,t,i)+ & gamz(1,1,k,t,i)*y(t-1,i)+gamz(2,1,k,t,i)*y(t-1,i)**2.d0+ & gamz(1,2,k,t,i)*y(t-2,i)+gamz(2,2,k,t,i)*y(t-2,i)**2.d0)/ & (1.d0+sumexp) 52 continue pc(K1)=1.d0/(1.d0+sumexp) c write(*,*)i,t,pc CALL RNMTN (1, 1, K1, real(pc), IR,1) do 53 k=1,K1 is=0 do ig1=1,k id(ig1,t,i)=ir(1,ig1) is=is+ir(1,ig1) enddo r(k,t,i)=is if(ir(1,k).eq.1) then y(t,i)=yv(k) y1(t,i)=dble(k)-1.d0 endif 53 continue endif 42 continue endif 40 continue inr=0 do i=1,N inr=inr+nr(i) c write(*,*)'i=',i c write(*,*)(y1(t,i),t=1,nr(i)) enddo write(*,*)inr c read(*,*)iii call findmle(x,y,yv,z,r,id,nr,b1,ba11, & N,IT,K1,nx,nalpha1,ialph,ifail) if(ifail.eq.0) then write(*,75)(b1(ig),ig=1,K1-1+nx-1) write(3,75)(b1(ig),ig=1,K1-1+nx-1) write(*,75)(ba11(ig),ig=1,(K1-1)*nalpha1*ialph) 75 format(100f12.4) endif 70 continue close(3) stop end c Subroutine for calculating \triangle_{itk} c for findcand subroutine caltri1(yv,z,pi,pj,ralpha0,gamz,N,nr,IT,K1, & ialph,nalpha,norder,triangle,kk) integer nr(N),t double precision ralpha0(nalpha,ialph,K1-1,norder), & triangle(K1-1,IT,N),tri(K1-1),tri0(K1-1), & z(nalpha,IT,N), & gamz(ialph,norder,K1-1,IT,N),yv(K1),temptri(K1-1), & f(K1-1),df(K1-1,K1-1),ainv(K1-1,K1-1), & pc(K1,K1,K1),pi(K1,IT,N),dpc(K1-1,K1,K1,K1-1), & pj(K1,K1,IT,N), & epsilon,sum0,temp0,tempsum epsilon=0.0001d0 no1=kk do 90 l=1,norder do 90 j=1,ialph do 92 i=1,N if(nr(i).ge.2) then do 93 t=2,nr(i) do 95 k=1,K1-1 sum0=0.d0 do 97 ig=1,nalpha sum0=sum0+z(ig,t,i)*ralpha0(ig,j,k,l) 97 continue gamz(j,l,k,t,i)=sum0 if((t.eq.2).and.(l.eq.2)) then gamz(j,l,k,t,i)=0.d0 endif 95 continue 93 continue endif 92 continue 90 continue c write(*,*)'In the caltri0' do 100 i=1,N if(nr(i).ge.2) then do 103 t=2,nr(i) if(kk.eq.1) then write(*,*)i,t endif if(no1.eq.1)then call DRNUN(K1-1,tri) endif do 105 ig=1,K1-1 if(kk.ge.2)tri(ig)=triangle(ig,t,i) if(kk.eq.1) then tri(ig)=tri(ig)-0.5d0 endif 105 continue niter=10000 if(t.eq.2) then do 107 no=1,niter do 108 ig=1,K1-1 tri0(ig)=tri(ig) 108 continue do 110 ik=1,K1-1 do 110 ig=1,K1-1 df(ik,ig)=0.d0 110 continue do 112 k=1,K1-1 f(k)=0.d0 do 114 j=1,K1 tempsum=0.d0 do 116 l=1,K1-1 tempsum=tempsum+dexp(tri(l)+gamz(1,1,l,t,i)*yv(j)+ & gamz(2,1,l,t,i)*yv(j)**2.d0) 116 continue do 118 ik=1,K1 if(ik.lt.K1) then pc(1,j,ik)=dexp(tri(ik)+gamz(1,1,ik,t,i)*yv(j)+ & gamz(2,1,ik,t,i)*yv(j)**2.d0)/(1.d0+tempsum) else pc(1,j,K1)=1.d0/(1.d0+tempsum) endif 118 continue do 120 ig=1,K1-1 dpc(ig,1,j,k)=0.d0 if(ig.eq.k) dpc(ig,1,j,k)=pc(1,j,k)*(1.d0-pc(1,j,k)) if(ig.ne.k) dpc(ig,1,j,k)=-pc(1,j,k)*pc(1,j,ig) 120 continue f(k)=f(k)+pc(1,j,k)*pi(j,t-1,i) do 122 ig=1,K1-1 df(k,ig)=df(k,ig)+dpc(ig,1,j,k)*pi(j,t-1,i) 122 continue 114 continue f(k)=f(k)-pi(k,t,i) 112 continue call DLINRG(K1-1,df,K1-1,ainv,K1-1) call DMURRV(K1-1,K1-1,ainv,K1-1,K1-1,f,1,K1-1,temptri) do 123 ig=1,K1-1 tri(ig)=tri0(ig)-temptri(ig)/2.d0 123 continue temp0=0.d0 do 124 ig=1,K1-1 temp0=temp0+abs(tri(ig)-tri0(ig)) 124 continue if(temp0.lt.epsilon) then do 126 j1=1,K1 do 126 j2=1,K1 pj(j2,j1,t,i)=pc(1,j2,j1)*pi(j2,t-1,i) 126 continue goto 150 endif do 125 ig=1,K1-1 if(abs(tri(ig)).gt.30) then call DRNUN(K1-1,tri) goto 107 endif 125 continue if(no.eq.niter) write(*,*)'Not converge in triangle' 107 continue else do 130 no=1,niter do 132 ig=1,K1-1 tri0(ig)=tri(ig) 132 continue do 134 ik=1,K1-1 do 134 ig=1,K1-1 df(ik,ig)=0.d0 134 continue do 136 k=1,K1-1 f(k)=0.d0 do 137 j1=1,K1 do 137 j2=1,K1 tempsum=0.d0 do 138 l=1,K1-1 tempsum=tempsum+dexp(tri(l)+gamz(1,1,l,t,i)*yv(j1)+ & gamz(2,1,l,t,i)*yv(j1)**2.d0+gamz(1,2,l,t,i)*yv(j2)+ & gamz(2,2,l,t,i)*yv(j2)**2.d0) 138 continue do 139 ik=1,K1 if(ik.lt.K1) then pc(j2,j1,ik)=dexp(tri(ik)+gamz(1,1,ik,t,i)*yv(j1)+ & gamz(2,1,ik,t,i)*yv(j1)**2.d0+gamz(1,2,ik,t,i)*yv(j2)+ & gamz(2,2,ik,t,i)*yv(j2)**2.d0)/(1.d0+tempsum) else pc(j2,j1,K1)=1.d0/(1.d0+tempsum) endif 139 continue do 140 ig=1,K1-1 dpc(ig,j2,j1,k)=0.d0 if(ig.eq.k) dpc(ig,j2,j1,k)= & pc(j2,j1,k)*(1.d0-pc(j2,j1,k)) if(ig.ne.k) dpc(ig,j2,j1,k)=-pc(j2,j1,k)*pc(j2,j1,ig) 140 continue f(k)=f(k)+pc(j2,j1,k)*pj(j2,j1,t-1,i) do 141 ig=1,K1-1 df(k,ig)=df(k,ig)+dpc(ig,j2,j1,k)*pj(j2,j1,t-1,i) 141 continue 137 continue f(k)=f(k)-pi(k,t,i) 136 continue call DLINRG(K1-1,df,K1-1,ainv,K1-1) call DMURRV(K1-1,K1-1,ainv,K1-1,K1-1,f,1,K1-1,temptri) do 142 ig=1,K1-1 tri(ig)=tri0(ig)-temptri(ig)/4.d0 142 continue temp0=0.d0 do 143 ig=1,K1-1 temp0=temp0+abs(tri(ig)-tri0(ig)) 143 continue if(temp0.lt.epsilon) then do 146 j1=1,K1 do 146 j2=1,K1 pj(j2,j1,t,i)=0.d0 do 147 ig=1,K1 pj(j2,j1,t,i)=pj(j2,j1,t,i)+pc(ig,j2,j1)*pj(ig,j2,t-1,i) 147 continue 146 continue goto 150 endif do 145 ig=1,K1-1 if(abs(tri(ig)).gt.30) then call DRNUN(K1-1,tri) goto 130 endif 145 continue if(no.eq.niter) write(*,*)'Not converge in triangle' 130 continue endif 150 no1=1 do 154 k=1,K1-1 triangle(k,t,i)=tri(k) 154 continue c write(*,*)(tri(k),k=1,K1-1) 103 continue endif 100 continue return end c Subroutine to calculate p^M subroutine calpm(beta1,x,N,nr,IT,K1,nx,pm) integer nr(N),t double precision beta1(nx,K1),x(nx,IT,N),pm(K1,IT,N),temp do 153 i=1,N do 153 t=1,nr(i) do 155 k=1,K1 temp=0.d0 do 157 ig=1,nx temp=temp+x(ig,t,i)*beta1(ig,k) 157 continue if(k.lt.K1) then pm(k,t,i)=dexp(temp)/(1.d0+dexp(temp)) else pm(K1,t,i)=1.d0 endif 155 continue c write(*,*)i,t,(pm(ik,t,i),ik=1,K1) c read(*,*)ia 153 continue return end C Subroutine to calculate pi subroutine calpi(pm,N,nr,IT,K1,pi) integer nr(N),t double precision pm(K1,IT,N),pi(K1,IT,N) do 170 i=1,N do 170 t=1,nr(i) do 175 k=1,K1 if(k.eq.1) then pi(k,t,i)=pm(k,t,i) else pi(k,t,i)=pm(k,t,i)-pm(k-1,t,i) endif 175 continue c write(*,*)(pi(ig,t,i),ig=1,K1) 170 continue return end c Subroutine for calculating \triangle_{itk} c for findcand subroutine caltri0(yv,z,pi,ralpha0,N,nr,IT,K1,ialph,nalpha, & triangle,kk) integer nr(N),t double precision ralpha0(nalpha,ialph,K1-1),triangle(K1-1,IT,N), & tri(K1-1),tri0(K1-1), & z(nalpha,IT,N), & gamz(ialph,K1-1,IT,N),yv(K1),temptri(K1-1), & f(K1-1),df(K1-1,K1-1),ainv(K1-1,K1-1), & pc(K1,K1),pi(K1,IT,N),dpc(K1-1,K1,K1-1), & epsilon,sum0 epsilon=0.001d0 no1=kk do 90 j=1,ialph do 92 i=1,N if(nr(i).ge.2) then do 93 t=2,nr(i) do 95 k=1,K1-1 sum0=0.d0 do 97 ig=1,nalpha sum0=sum0+z(ig,t,i)*ralpha0(ig,j,k) 97 continue gamz(j,k,t,i)=sum0 95 continue 93 continue endif 92 continue 90 continue do 100 i=1,N if(nr(i).ge.2) then do 103 t=2,nr(i) c if(kk.eq.1) then c write(*,*)i,t c endif if(no1.eq.1)then call DRNUN(K1-1,tri) endif do 105 ig=1,K1-1 if(kk.ge.2)tri(ig)=triangle(ig,t,i) if(kk.eq.1) then tri(ig)=tri(ig)-0.5d0 endif 105 continue niter=20000 do 107 no=1,niter do 110 ig=1,K1-1 tri0(ig)=tri(ig) 110 continue do 115 ik=1,K1-1 do 115 ig=1,K1-1 115 df(ik,ig)=0.d0 do 118 k=1,K1-1 f(k)=0.d0 do 120 j=1,K1 tempsum=0.d0 do 125 l=1,K1-1 tempsum=tempsum+dexp(tri(l)+gamz(1,l,t,i)*yv(j)+ & gamz(2,l,t,i)*yv(j)**2.d0) 125 continue do 128 ik=1,K1 if(ik.lt.K1) then pc(j,ik)=dexp(tri(ik)+gamz(1,ik,t,i)*yv(j)+ & gamz(2,ik,t,i)*yv(j)**2.d0)/(1.d0+tempsum) else pc(j,K1)=1.d0/(1.d0+tempsum) endif 128 continue do 130 ig=1,K1-1 if(ig.eq.k) dpc(ig,j,k)=pc(j,k)*(1-pc(j,k)) if(ig.ne.k) dpc(ig,j,k)=-pc(j,k)*pc(j,ig) 130 continue f(k)=f(k)+pc(j,k)*pi(j,t-1,i) do 133 ig=1,K1-1 df(k,ig)=df(k,ig)+dpc(ig,j,k)*pi(j,t-1,i) 133 continue 120 continue f(k)=f(k)-pi(k,t,i) 118 continue c do ig=1,(K1-1) c write(*,*)(df(ig,ik),ik=1,(K1-1)) c enddo call DLINRG(K1-1,df,K1-1,ainv,K1-1) call DMURRV(K1-1,K1-1,ainv,K1-1,K1-1,f,1,K1-1,temptri) do 135 ig=1,K1-1 tri(ig)=tri0(ig)-temptri(ig)/4.d0 135 continue temp0=0.d0 do 140 ig=1,K1-1 temp0=temp0+abs(tri(ig)-tri0(ig)) 140 continue if(temp0.lt.epsilon) goto 147 do 145 ig=1,K1-1 if(abs(tri(ig)).gt.30) then call DRNUN(K1-1,tri) goto 107 endif 145 continue if(no.eq.niter) write(*,*)'Not converge in triangle' 107 continue 147 no1=1 do 150 k=1,K1-1 triangle(k,t,i)=tri(k) 150 continue 103 continue endif 100 continue return end c----------------------------------------------------------c c Find Mle for beta_0, beta, alpha c c----------------------------------------------------------c subroutine findmle(x,y,yv,z,r,id,nr,b1,ba1, & N,IT,K1,nx,nalpha,ialph,ifail) integer nr(N),t,r(K1,IT,N),id(K1,IT,N) double precision x(nx,IT,N),z(nalpha,IT,N),yv(K1),y(IT,N), & b0(K1-1+nx-1),b1(K1-1+nx-1), & ba0((K1-1)*nalpha*ialph),ba1((K1-1)*nalpha*ialph), & ralpha0(nalpha,ialph,K1-1), & triangle(K1-1,IT,N),pm(K1,IT,N),pi(K1,IT,N), & beta1(nx,K1),alpsum1(K1-1),alpsum2(K1-1), & dpmdbet0(K1-1,K1-1,IT,N), dpmdbet(nx-1,K1-1,IT,N), & dtrdbet0(K1-1,K1-1,IT,N), dtrdbet(nx-1,K1-1,IT,N), & dtrdalp1(nalpha,K1-1,K1-1,IT,N), & dtrdalp2(nalpha,K1-1,K1-1,IT,N), & dhdgam1(K1-1,K1,K1,IT,N),dhdgam2(K1-1,K1,K1,IT,N), & sumalp1(nalpha,K1-1),sumalp2(nalpha,K1-1), & suma1(nalpha,nalpha,K1-1,K1-1),suma2(nalpha,nalpha,K1-1,K1-1), & suma21(nalpha,nalpha,K1-1,K1-1), & dQ1(K1-1+nx-1+(K1-1)*nalpha*ialph), & H1(K1-1+nx-1+(K1-1)*nalpha*ialph, & K1-1+nx-1+(K1-1)*nalpha*ialph), & H1inv(K1-1+nx-1+(K1-1)*nalpha*ialph, & K1-1+nx-1+(K1-1)*nalpha*ialph), & pc(K1,IT,N),h(K1,K1,IT,N),dh(K1-1,K1,K1,IT,N), & tempb1(K1-1+nx-1+(K1-1)*nalpha*ialph), & dl(K1-1+nx-1),ddb1(nx-1,nx-1), & ddb0(K1-1,K1-1),ddb10(nx-1,K1-1), & sumexph(K1),sumbet(nx-1),sumbet0(K1-1), & sumb1(nx-1,nx-1),sumb0(K1-1,K1-1),sumb10(nx-1,K1-1), & suma1b1(nalpha,nx-1,K1-1),suma1b0(nalpha,K1-1,K1-1), & suma2b1(nalpha,nx-1,K1-1),suma2b0(nalpha,K1-1,K1-1), & aloglik,aloglik0,tempdh,temp0,temp1,temp2,update do 480 ig=1,K1-1+nx-1 b0(ig)=b1(ig) 480 continue do 481 ig=1,(K1-1)*nalpha*ialph ba0(ig)=ba1(ig) 481 continue do 482 j=1,ialph do 482 k=1,K1-1 do 482 ig=1,nalpha ralpha0(ig,j,k)= & ba1((K1-1)*nalpha*(j-1)+(k-1)*nalpha+ig) 482 continue do 515 k=1,K1 if(k.lt.K1) then beta1(1,k)=b1(k) do 520 ig=1,nx-1 beta1(ig+1,k)=b1(K1-1+ig) 520 continue else do 525 ig=1,nx beta1(ig,K1)=0.d0 525 continue endif 515 continue nolog=0 do 508 iter=1,1000 c write(*,*)'Iteration=',iter c Calculate P^M call calpm(beta1,x,N,nr,IT,K1,nx,pm) c Calculate pi call calpi(pm,N,nr,IT,K1,pi) c Calculate \triangle_{itk} given \beta_0^(n-1), \beta^(n-1) call caltri0(yv,z,pi,ralpha0,N,nr,IT,K1,ialph,nalpha, & triangle,iter) aloglik0=aloglik aloglik=0.d0 c Calculate P_{itk}^c do 527 i=1,N do 526 ii=1,K1 aloglik=aloglik+dble(id(ii,1,i))*dlog(pi(ii,1,i)) 526 continue if(nr(i).ge.2) then do 530 t=2,nr(i) do 533 ik=1,K1-1 alpsum1(ik)=0.d0 alpsum2(ik)=0.d0 do 531 ig=1,nalpha alpsum1(ik)=alpsum1(ik)+z(ig,t,i)*ralpha0(ig,1,ik) alpsum2(ik)=alpsum2(ik)+z(ig,t,i)*ralpha0(ig,2,ik) 531 continue 533 continue sumexp=0.d0 do 535 l=1,K1-1 sumexp=sumexp+dexp(triangle(l,t,i)+alpsum1(l)*y(t-1,i)+ & alpsum2(l)*y(t-1,i)**2.d0) 535 continue do 537 k=1,K1-1 pc(k,t,i)=dexp(triangle(k,t,i)+alpsum1(k)*y(t-1,i)+ & alpsum2(k)*y(t-1,i)**2.d0)/(1.d0+sumexp) 537 continue pc(K1,t,i)=1.d0/(1.d0+sumexp) do 534 ik1=1,K1-1 aloglik=aloglik+dble(id(ik1,t,i))* & (triangle(ik1,t,i)+alpsum1(ik1)*y(t-1,i)+ & alpsum2(ik1)*y(t-1,i)**2.d0) 534 continue aloglik=aloglik+dlog(pc(K1,t,i)) do 538 j=1,K1 sumexph(j)=0.d0 do 540 l=1,K1-1 sumexph(j)=sumexph(j)+dexp(triangle(l,t,i)+ & alpsum1(l)*yv(j)+alpsum2(l)*yv(j)**2.d0) 540 continue do 543 k=1,K1-1 h(j,k,t,i)=dexp(triangle(k,t,i)+alpsum1(k)*yv(j)+ & alpsum2(k)*yv(j)**2.d0)/(1.d0+sumexph(j)) 543 continue h(j,K1,t,i)=1.d0/(1.d0+sumexph(j)) 538 continue 530 continue endif 527 continue write(*,*)'loglikelihood =',aloglik if(aloglik.lt.aloglik0) then nolog=nolog+1 if(iter.ge.3) then if((nolog.gt.35).or.(dabs(aloglik-aloglik0).gt.50)) then ifail=1 goto 860 endif endif endif c Calculate d P^M / d beta_0 and d P^M / d beta do 545 i=1,N do 545 t=1,nr(i) do 545 k=1,K1-1 do 546 j=1,K1-1 if(k.eq.j) dpmdbet0(j,k,t,i)=pm(k,t,i)*(1-pm(k,t,i)) if(k.ne.j) dpmdbet0(j,k,t,i)=0.d0 546 continue do 547 ig=1,nx-1 dpmdbet(ig,k,t,i)=x(ig+1,t,i)*pm(k,t,i)*(1.d0-pm(k,t,i)) 547 continue 545 continue c Calculate d h_{ikg}^(t) / d gamma_{it1}^(a) and d h_{ikg}^(t) / d gamma_{it2}^(b) do 548 i=1,N if(nr(i).ge.2) then do 549 t=2,nr(i) do 549 k=1,K1 do 549 ig=1,K1 do 549 ia=1,K1-1 if(k.eq.ia)then dhdgam1(ia,ig,k,t,i)=h(ig,k,t,i)*(1-h(ig,k,t,i))*yv(ig) dhdgam2(ia,ig,k,t,i)= & h(ig,k,t,i)*(1-h(ig,k,t,i))*yv(ig)**2.d0 endif if(k.ne.ia)then dhdgam1(ia,ig,k,t,i)=-h(ig,k,t,i)*h(ig,ia,t,i)*yv(ig) dhdgam2(ia,ig,k,t,i)= & -h(ig,k,t,i)*h(ig,ia,t,i)*yv(ig)**2.d0 endif 549 continue endif 548 continue c Calculate d tri / d beta call dtrbet(z,ralpha0,triangle,yv,pi,dpmdbet, & N,IT,nr,K1,nalpha,nx,ialph,dtrdbet) c Calculate d tri / d beta_0 call dtrbet0(z,ralpha0,triangle,yv,pi,dpmdbet0,N,IT, & nr,K1,nalpha,ialph,dtrdbet0) c Calculate d tri / d alpha_1^(a) call dtralp1(z,ralpha0,triangle,yv,pi,dhdgam1,N,IT, & nr,K1,nalpha,ialph,dtrdalp1) c Calculate d tri / d alpha_2^(b) call dtralp2(z,ralpha0,triangle,yv,pi,dhdgam2,N,IT, & nr,K1,nalpha,ialph,dtrdalp2) c Calculate Information matrix at time=1 call findhess1(r,N,IT,K1,nx,pm,dpmdbet0,dpmdbet,dl, & ddb1,ddb0,ddb10) do 550 ii=1,nx-1 sumbet(ii)=0.d0 550 continue do 551 j=1,K1-1 sumbet0(j)=0.d0 551 continue do 552 ia=1,K1-1 do 552 j=1,nalpha sumalp1(j,ia)=0.d0 sumalp2(j,ia)=0.d0 552 continue do 553 j=1,K1-1 do 553 i=1,N if(nr(i).ge.2) then do 554 t=2,nr(i) do 555 k=1,K1-1 sumbet0(j)=sumbet0(j)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdbet0(j,k,t,i) 555 continue 554 continue endif 553 continue do 557 ia=1,nx-1 do 557 i=1,N if(nr(i).ge.2) then do 559 t=2,nr(i) do 560 k=1,K1-1 sumbet(ia)=sumbet(ia)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdbet(ia,k,t,i) 560 continue 559 continue endif 557 continue do 563 ia=1,K1-1 do 563 ig=1,nalpha do 564 i=1,N if(nr(i).ge.2) then do 565 t=2,nr(i) do 567 k=1,K1-1 sumalp1(ig,ia)=sumalp1(ig,ia)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdalp1(ig,ia,k,t,i) sumalp2(ig,ia)=sumalp2(ig,ia)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdalp2(ig,ia,k,t,i) 567 continue sumalp1(ig,ia)=sumalp1(ig,ia)+ & (dble(id(ia,t,i))-pc(ia,t,i))*y(t-1,i)*z(ig,t,i) sumalp2(ig,ia)=sumalp2(ig,ia)+ & (dble(id(ia,t,i))-pc(ia,t,i))*y(t-1,i)**2.d0*z(ig,t,i) 565 continue endif 564 continue 563 continue do 568 ig=1,K1-1 dQ1(ig)=sumbet0(ig)+dl(ig) 568 continue do 569 ig=1,nx-1 dQ1(ig+K1-1)=sumbet(ig)+dl(K1-1+ig) 569 continue do 570 ia=1,K1-1 do 570 ig=1,nalpha dQ1(K1-1+nx-1+(ia-1)*nalpha+ig)=sumalp1(ig,ia) dQ1(K1-1+nx-1+(K1-1)*nalpha+(ia-1)*nalpha+ig)=sumalp2(ig,ia) 570 continue c Calculate Hessian matrix do 571 ia1=1,K1-1 do 571 ia2=1,K1-1 sumb0(ia1,ia2)=0.d0 571 continue do 573 ig1=1,nx-1 do 573 ig2=1,nx-1 sumb1(ig1,ig2)=0.d0 573 continue do 574 ig=1,nx-1 do 574 ia=1,K1-1 sumb10(ig,ia)=0.d0 574 continue do 575 i1=1,K1-1 do 575 i2=1,K1-1 do 575 ig1=1,nalpha do 575 ig2=1,nalpha suma1(ig1,ig2,i1,i2)=0.d0 suma2(ig1,ig2,i1,i2)=0.d0 suma21(ig1,ig2,i1,i2)=0.d0 575 continue do 580 i=1,N if(nr(i).ge.2)then do 581 t=2,nr(i) do 581 k=1,K1 do 583 ig=1,K1-1 do 585 j=1,K1 if(ig.eq.k) then dh(ig,j,k,t,i)=h(j,k,t,i)*(1.d0-h(j,k,t,i)) else if(ig.ne.k) then dh(ig,j,k,t,i)=-h(j,k,t,i)*h(j,ig,t,i) endif 585 continue 583 continue 581 continue endif 580 continue c hessian matrix for \beta_0j \beta_0j' do 587 i1=1,K1-1 do 587 i2=1,K1-1 do 590 i=1,N if(nr(i).ge.2) then do 591 t=2,nr(i) do 593 k=1,K1-1 do 594 ig=1,K1-1 tempdh=0.d0 do 595 ik=1,K1 tempdh=tempdh+dh(ig,ik,k,t,i)*pi(ik,t-1,i) 595 continue sumb0(i1,i2)=sumb0(i1,i2)+ & tempdh*dtrdbet0(i1,ig,t,i)*dtrdbet0(i2,k,t,i) 594 continue 593 continue 591 continue endif 590 continue 587 continue c \beta \beta do 600 i1=1,nx-1 do 600 i2=1,nx-1 do 601 i=1,N if(nr(i).ge.2) then do 603 t=2,nr(i) do 605 k=1,K1-1 do 607 ig=1,K1-1 tempdh=0.d0 do 608 ik=1,K1 tempdh=tempdh+dh(ig,ik,k,t,i)*pi(ik,t-1,i) 608 continue sumb1(i1,i2)=sumb1(i1,i2)+ & tempdh*dtrdbet(i1,ig,t,i)*dtrdbet(i2,k,t,i) 607 continue 605 continue 603 continue endif 601 continue 600 continue c \beta \beta_0j do 610 i1=1,nx-1 do 610 i2=1,K1-1 do 611 i=1,N if(nr(i).ge.2) then do 613 t=2,nr(i) do 615 k=1,K1-1 do 617 ig=1,K1-1 tempdh=0.d0 do 620 ik=1,K1 tempdh=tempdh+dh(ig,ik,k,t,i)*pi(ik,t-1,i) 620 continue sumb10(i1,i2)=sumb10(i1,i2)+ & tempdh*dtrdbet(i1,ig,t,i)*dtrdbet0(i2,k,t,i) 617 continue 615 continue 613 continue endif 611 continue c write(*,*)'beta,beta_0=',sumb10(i1,i2),i1,i2 610 continue c \alpha_1 \alph_1 do 623 i1=1,K1-1 do 623 i2=1,K1-1 do 623 ig1=1,nalpha do 623 ig2=1,nalpha do 625 i=1,N if(nr(i).ge.2) then do 626 t=2,nr(i) do 628 k=1,K1-1 do 630 ig=1,K1-1 tempdh=0.d0 do 632 ij=1,K1 tempdh=tempdh+dh(ig,ij,k,t,i)*pi(ij,t-1,i) 632 continue suma1(ig1,ig2,i1,i2)=suma1(ig1,ig2,i1,i2)+ & tempdh*dtrdalp1(ig1,i1,ig,t,i)*dtrdalp1(ig2,i2,k,t,i) 630 continue temp0=0.d0 do 635 ij=1,K1 temp0=temp0+dhdgam1(i1,ij,k,t,i)*pi(ij,t-1,i) 635 continue suma1(ig1,ig2,i1,i2)=suma1(ig1,ig2,i1,i2)+ & temp0*z(ig1,t,i)*dtrdalp1(ig2,i2,k,t,i) 628 continue do 637 ig=1,K1-1 temp1=0.d0 do 640 ij=1,K1 temp1=temp1+dh(ig,ij,i2,t,i)*yv(ij)*pi(ij,t-1,i) 640 continue suma1(ig1,ig2,i1,i2)=suma1(ig1,ig2,i1,i2)+ & temp1*dtrdalp1(ig1,i1,ig,t,i)*z(ig2,t,i) 637 continue temp2=0.d0 do 643 ij=1,K1 temp2=temp2+dhdgam1(i1,ij,i2,t,i)*yv(ij)*pi(ij,t-1,i) 643 continue suma1(ig1,ig2,i1,i2)=suma1(ig1,ig2,i1,i2)+ & temp2*z(ig1,t,i)*z(ig2,t,i) 626 continue endif 625 continue 623 continue c \alpha_2 \alph_2 do 650 i1=1,K1-1 do 650 i2=1,K1-1 do 650 ig1=1,nalpha do 650 ig2=1,nalpha do 653 i=1,N if(nr(i).ge.2) then do 655 t=2,nr(i) do 657 k=1,K1-1 do 658 ig=1,K1-1 tempdh=0.d0 do 660 ij=1,K1 tempdh=tempdh+dh(ig,ij,k,t,i)*pi(ij,t-1,i) 660 continue suma2(ig1,ig2,i1,i2)=suma2(ig1,ig2,i1,i2)+ & tempdh*dtrdalp2(ig1,i1,ig,t,i)*dtrdalp2(ig2,i2,k,t,i) 658 continue temp0=0.d0 do 663 ij=1,K1 temp0=temp0+dhdgam2(i1,ij,k,t,i)*pi(ij,t-1,i) 663 continue suma2(ig1,ig2,i1,i2)=suma2(ig1,ig2,i1,i2)+ & temp0*z(ig1,t,i)*dtrdalp2(ig2,i2,k,t,i) 657 continue do 665 ig=1,K1-1 temp1=0.d0 do 667 ij=1,K1 temp1=temp1+dh(ig,ij,i2,t,i)*yv(ij)**2*pi(ij,t-1,i) 667 continue suma2(ig1,ig2,i1,i2)=suma2(ig1,ig2,i1,i2)+ & temp1*dtrdalp2(ig1,i1,ig,t,i)*z(ig2,t,i) 665 continue temp2=0.d0 do 670 ij=1,K1 temp2=temp2+dhdgam2(i1,ij,i2,t,i)*yv(ij)**2*pi(ij,t-1,i) 670 continue suma2(ig1,ig2,i1,i2)=suma2(ig1,ig2,i1,i2)+ & temp2*z(ig1,t,i)*z(ig2,t,i) c write(*,*)'sum2=',suma2(ig1,ig2,i1,i2),i1,i2 655 continue endif 653 continue 650 continue c \alpha_2 \alph_1 do 680 i1=1,K1-1 do 680 i2=1,K1-1 do 680 ig1=1,nalpha do 680 ig2=1,nalpha do 683 i=1,N if(nr(i).ge.2) then do 685 t=2,nr(i) do 687 k=1,K1-1 do 688 ig=1,K1-1 tempdh=0.d0 do 690 ij=1,K1 tempdh=tempdh+dh(ig,ij,k,t,i)*pi(ij,t-1,i) 690 continue suma21(ig1,ig2,i1,i2)=suma21(ig1,ig2,i1,i2)+ & tempdh*dtrdalp2(ig1,i1,ig,t,i)*dtrdalp1(ig2,i2,k,t,i) 688 continue temp0=0.d0 do 691 ij=1,K1 temp0=temp0+dhdgam2(i1,ij,k,t,i)*pi(ij,t-1,i) 691 continue suma21(ig1,ig2,i1,i2)=suma21(ig1,ig2,i1,i2)+ & temp0*z(ig1,t,i)*dtrdalp1(ig2,i2,k,t,i) 687 continue do 693 ig=1,K1-1 temp1=0.d0 do 695 ij=1,K1 temp1=temp1+dh(ig,ij,i2,t,i)*yv(ij)*pi(ij,t-1,i) 695 continue suma21(ig1,ig2,i1,i2)=suma21(ig1,ig2,i1,i2)+ & temp1*dtrdalp2(ig1,i1,ig,t,i)*z(ig2,t,i) 693 continue temp2=0.d0 do 697 ij=1,K1 temp2=temp2+dhdgam2(i1,ij,i2,t,i)*yv(ij)*pi(ij,t-1,i) 697 continue suma21(ig1,ig2,i1,i2)=suma21(ig1,ig2,i1,i2)+ & temp2*z(ig1,t,i)*z(ig2,t,i) 685 continue endif 683 continue 680 continue c alpha_1 beta and alpha_2 beta do 700 i2=1,K1-1 do 700 ig1=1,nx-1 do 700 ig2=1,nalpha do 702 i=1,N if(nr(i).ge.2) then do 703 t=2,nr(i) do 704 k=1,K1-1 do 705 ig=1,K1-1 tempdh=0.d0 do 706 ik=1,K1 tempdh=tempdh+dh(ig,ik,k,t,i)*pi(ik,t-1,i) 706 continue suma1b1(ig2,ig1,i2)=suma1b1(ig2,ig1,i2)+ & tempdh*dtrdalp1(ig2,i2,ig,t,i)*dtrdbet(ig1,k,t,i) suma2b1(ig2,ig1,i2)=suma2b1(ig2,ig1,i2)+ & tempdh*dtrdalp2(ig2,i2,ig,t,i)*dtrdbet(ig1,k,t,i) 705 continue temp1=0.d0 do 710 ij=1,K1 temp1=temp1+dhdgam1(i2,ij,k,t,i)*pi(ij,t-1,i) 710 continue suma1b1(ig2,ig1,i2)=suma1b1(ig2,ig1,i2)+ & temp1*z(ig2,t,i)*dtrdbet(ig1,k,t,i) temp2=0.d0 do 711 ij=1,K1 temp2=temp2+dhdgam2(i2,ij,k,t,i)*pi(ij,t-1,i) 711 continue suma2b1(ig2,ig1,i2)=suma2b1(ig2,ig1,i2)+ & temp2*z(ig2,t,i)*dtrdbet(ig1,k,t,i) 704 continue 703 continue endif 702 continue 700 continue c \alpha_1 \beta_0 and \alpha_2 \beta_0 do 712 i2=1,K1-1 do 712 ig1=1,K1-1 do 712 ig2=1,nalpha do 713 i=1,N if(nr(i).ge.2) then do 714 t=2,nr(i) do 715 k=1,K1-1 do 716 ig=1,K1-1 tempdh=0.d0 do 718 ik=1,K1 tempdh=tempdh+dh(ig,ik,k,t,i)*pi(ik,t-1,i) 718 continue suma1b0(ig2,ig1,i2)=suma1b0(ig2,ig1,i2)+ & tempdh*dtrdalp1(ig2,i2,ig,t,i)*dtrdbet0(ig1,k,t,i) suma2b0(ig2,ig1,i2)=suma2b0(ig2,ig1,i2)+ & tempdh*dtrdalp2(ig2,i2,ig,t,i)*dtrdbet0(ig1,k,t,i) 716 continue temp1=0.d0 do 720 ij=1,K1 temp1=temp1+dhdgam1(i2,ij,k,t,i)*pi(ij,t-1,i) 720 continue suma1b0(ig2,ig1,i2)=suma1b0(ig2,ig1,i2)+ & temp1*z(ig2,t,i)*dtrdbet0(ig1,k,t,i) temp2=0.d0 do 722 ij=1,K1 temp2=temp2+dhdgam2(i2,ij,k,t,i)*pi(ij,t-1,i) 722 continue suma2b0(ig2,ig1,i2)=suma2b0(ig2,ig1,i2)+ & temp2*z(ig2,t,i)*dtrdbet0(ig1,k,t,i) 715 continue 714 continue endif 713 continue 712 continue c construct Hessain matrix do 800 ik=1,K1-1 do 803 ig=1,K1-1 H1(ik,ig)=sumb0(ik,ig)+ddb0(ik,ig) 803 continue do 805 ig=1,nx-1 H1(ig+K1-1,ik)=sumb10(ig,ik)+ddb10(ig,ik) H1(ik,K1-1+ig)=sumb10(ig,ik)+ddb10(ig,ik) 805 continue 800 continue do 807 ik=1,nx-1 do 807 ig=1,nx-1 H1(K1-1+ik,K1-1+ig)=sumb1(ik,ig)+ddb1(ik,ig) 807 continue do 810 ik1=1,K1-1 do 810 ik2=1,K1-1 do 812 ig1=1,nalpha do 812 ig2=1,nalpha H1(K1-1+nx-1+(ik1-1)*nalpha+ig1, & K1-1+nx-1+(ik2-1)*nalpha+ig2)=suma1(ig1,ig2,ik1,ik2) H1(K1-1+nx-1+(K1-1)*nalpha+(ik1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*nalpha+(ik2-1)*nalpha+ig2)= & suma2(ig1,ig2,ik1,ik2) 812 continue do 814 ig1=1,nalpha do 814 ig2=1,nalpha H1(K1-1+nx-1+(K1-1)*nalpha+(ik1-1)*nalpha+ig1, & K1-1+nx-1+(ik2-1)*nalpha+ig2)=suma21(ig1,ig2,ik1,ik2) H1(K1-1+nx-1+(ik2-1)*nalpha+ig2, & K1-1+nx-1+(K1-1)*nalpha+(ik1-1)*nalpha+ig1)= & suma21(ig1,ig2,ik1,ik2) 814 continue 810 continue do 816 ik1=1,K1-1 do 816 ig1=1,nalpha do 816 ig2=1,K1-1 H1(K1-1+nx-1+(ik1-1)*nalpha+ig1,ig2)=suma1b0(ig1,ig2,ik1) H1(ig2,K1-1+nx-1+(ik1-1)*nalpha+ig1)=suma1b0(ig1,ig2,ik1) 816 continue do 818 ik1=1,K1-1 do 818 ig1=1,nalpha do 818 ig2=1,K1-1 H1(K1-1+nx-1+(K1-1)*nalpha+(ik1-1)*nalpha+ig1,ig2)= & suma2b0(ig1,ig2,ik1) H1(ig2,K1-1+nx-1+(K1-1)*nalpha+(ik1-1)*nalpha+ig1)= & suma2b0(ig1,ig2,ik1) 818 continue do 820 ik2=1,K1-1 do 820 ig1=1,nx-1 do 820 ig2=1,nalpha H1(K1-1+nx-1+(ik2-1)*nalpha+ig2,K1-1+ig1)=suma1b1(ig2,ig1,ik2) H1(K1-1+ig1,K1-1+nx-1+(ik2-1)*nalpha+ig2)=suma1b1(ig2,ig1,ik2) 820 continue do 822 ik2=1,K1-1 do 822 ig1=1,nx-1 do 822 ig2=1,nalpha H1(K1-1+nx-1+(K1-1)*nalpha+(ik2-1)*nalpha+ig2,K1-1+ig1)= & suma2b1(ig2,ig1,ik2) H1(K1-1+ig1,K1-1+nx-1+(K1-1)*nalpha+(ik2-1)*nalpha+ig2)= & suma2b1(ig2,ig1,ik2) 822 continue call DLINRG(K1-1+nx-1+(K1-1)*nalpha*2,H1, & K1-1+nx-1+(K1-1)*nalpha*2,H1inv,K1-1+nx-1+(K1-1)*nalpha*2) call DMURRV(K1-1+nx-1+(K1-1)*nalpha*2, & K1-1+nx-1+(K1-1)*nalpha*2,H1inv,K1-1+nx-1+(K1-1)*nalpha*2, & K1-1+nx-1+(K1-1)*nalpha*2,dQ1,1,K1-1+nx-1+(K1-1)*nalpha*2, & tempb1) do 833 ig=1,K1-1+nx-1 b1(ig)=b0(ig)+tempb1(ig)/2.d0 833 continue do 834 ig=1,(K1-1)*nalpha*ialph ba1(ig)=ba0(ig)+tempb1(K1-1+nx-1+ig)/2.d0 834 continue update=0.d0 do 835 ig=1,K1-1+nx-1 update=update+abs(b1(ig)-b0(ig)) 835 continue do 836 ig=1,(K1-1)*nalpha*ialph update=update+abs(ba1(ig)-ba0(ig)) 836 continue if(update.le.0.0001d0) then ifail=0 goto 850 endif do 837 ig=1,K1-1+nx-1 b0(ig)=b1(ig) 837 continue do 838 ig=1,(K1-1)*nalpha*ialph ba0(ig)=ba1(ig) 838 continue do 840 j=1,ialph do 840 k=1,K1-1 do 840 ig=1,nalpha ralpha0(ig,j,k)= & ba1((K1-1)*nalpha*(j-1)+(k-1)*nalpha+ig) 840 continue do 843 k=1,K1 if(k.lt.K1) then beta1(1,k)=b1(k) do 845 ig=1,nx-1 beta1(ig+1,k)=b1(K1-1+ig) 845 continue else do 847 ig=1,nx beta1(ig,K1)=0.d0 847 continue endif 843 continue 508 continue 850 write(*,*)'Great! It is convergent' 860 return end c Calculate Information matrix at t=1 subroutine findhess1(r,N,IT,K1,nx,pm,dpmdbet0,dpmdbet,dl, & ddb1,ddb0,ddb10) c Reference : McCullagh(JRSS-B,1988) c Regression Models for Ordinal Data integer r(K1,IT,N) double precision pm(K1,IT,N), & dpmdbet0(K1-1,K1-1,IT,N), dpmdbet(nx-1,K1-1,IT,N), & dlike0(K1-1),dlike1(nx-1),dl(K1-1+nx-1), & ddb1(nx-1,nx-1), & ddb0(K1-1,K1-1),ddb10(nx-1,K1-1) c Calculate dlogL^(1)/d\beta and dlogL^(1)/d\beta_{0j} do 2060 j=1,K1-1 dlike0(j)=0.d0 2060 continue do 2063 j=1,K1-1 do 2063 i=1,N do 2063 k=1,K1-1 if(k.lt.K1-1) then dlike0(j)=dlike0(j)+(dble(r(k,1,i))-dble(r(k+1,1,i))* & pm(k,1,i)/pm(k+1,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet0(j,k,1,i)- & 1.d0/(pm(k+1,1,i)-pm(k,1,i))*dpmdbet0(j,k+1,1,i)) else dlike0(j)=dlike0(j)+(dble(r(k,1,i))-dble(r(k+1,1,i))* & pm(k,1,i)/pm(k+1,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet0(j,k,1,i)) endif 2063 continue do 2065 ig=1,nx-1 dlike1(ig)=0.d0 2065 continue do 2067 i=1,N do 2067 k=1,K1-1 do 2067 ig=1,nx-1 if(k.lt.K1-1) then dlike1(ig)=dlike1(ig)+ & (dble(r(k,1,i))-dble(r(k+1,1,i))*pm(k,1,i)/pm(k+1,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet(ig,k,1,i)- & 1.d0/(pm(k+1,1,i)-pm(k,1,i))*dpmdbet(ig,k+1,1,i)) else dlike1(ig)=dlike1(ig)+ & (dble(r(k,1,i))-dble(r(k+1,1,i))*pm(k,1,i)/pm(k+1,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet(ig,k,1,i)) endif 2067 continue do 2070 ik=1,nx-1 do 2070 ig=1,nx-1 ddb1(ik,ig)=0.d0 2070 continue do 2073 ik=1,nx-1 do 2073 ig=1,nx-1 do 2073 i=1,N do 2073 k=1,K1-1 if(k.lt.K1-1) then ddb1(ik,ig)=ddb1(ik,ig)+(dpmdbet(ik,k,1,i)-pm(k,1,i)/ & pm(k+1,1,i)*dpmdbet(ik,k+1,1,i))*(pm(k+1,1,i)/ & (pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet(ig,k,1,i)-1.d0/(pm(k+1,1,i)-pm(k,1,i))* & dpmdbet(ig,k+1,1,i)) else ddb1(ik,ig)=ddb1(ik,ig)+ & (dpmdbet(ik,k,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet(ig,k,1,i)) endif 2073 continue c write(*,*)'c4' do 2080 j1=1,K1-1 do 2080 j2=1,K1-1 ddb0(j1,j2)=0.d0 2080 continue do 2083 j1=1,K1-1 do 2083 j2=1,K1-1 do 2083 i=1,N do 2083 k=1,K1-1 if(k.lt.K1-1) then ddb0(j1,j2)=ddb0(j1,j2)+ & (dpmdbet0(j1,k,1,i)-pm(k,1,i)/pm(k+1,1,i)* & dpmdbet0(j1,k+1,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet0(j2,k,1,i)- & 1.d0/(pm(k+1,1,i)-pm(k,1,i))*dpmdbet0(j2,k+1,1,i)) else ddb0(j1,j2)=ddb0(j1,j2)+ & (dpmdbet0(j1,k,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet0(j2,k,1,i)) endif 2083 continue c write(*,*)ddbeta0 do 2085 j=1,K1-1 do 2085 ig=1,nx-1 ddb10(ig,j)=0.d0 2085 continue do 2087 j=1,K1-1 do 2087 ig=1,nx-1 do 2087 i=1,N do 2087 k=1,K1-1 if(k.lt.K1-1) then ddb10(ig,j)=ddb10(ig,j)+ & (dpmdbet(ig,k,1,i)-pm(k,1,i)/pm(k+1,1,i)* & dpmdbet(ig,k+1,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet0(j,k,1,i)- & 1.d0/(pm(k+1,1,i)-pm(k,1,i))*dpmdbet0(j,k+1,1,i)) else ddb10(ig,j)=ddb10(ig,j)+ & (dpmdbet(ig,k,1,i))* & (pm(k+1,1,i)/(pm(k,1,i)*(pm(k+1,1,i)-pm(k,1,i)))* & dpmdbet0(j,k,1,i)) endif 2087 continue c write(*,*)ddbeta10 do 2090 ig=1,K1-1 2090 dl(ig)=dlike0(ig) do 2093 ig=K1,nx-1+K1-1 2093 dl(ig)=dlike1(ig-(K1-1)) return end c Calculate d tri / d beta subroutine dtrbet(z,ralpha0,triangle,yv,pi,dpmdbet, & N,IT,nr,K1,nalpha,nx,ialph,dtrdbet) integer t,nr(N) double precision z(nalpha,IT,N), & ralpha0(nalpha,ialph,K1-1), & h0(K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,K1-1,IT,N), & triangle(K1-1,IT,N),yv(K1),dpmdbet(nx-1,K1-1,IT,N), & dhtsum(K1-1,K1-1),dtrdbet(nx-1,K1-1,IT,N),psum(nx-1,K1-1), & Ainv((K1-1)*(nx-1),(K1-1)*(nx-1)),res((K1-1)*(nx-1)), & dhtsum1((nx-1)*(K1-1),(nx-1)*(K1-1)), & psum1((K1-1)*(nx-1)),sum0,temp0,dht0 do 2100 j=1,ialph do 2102 i=1,N if(nr(i).ge.2) then do 2103 t=2,nr(i) do 2104 k=1,K1-1 sum0=0.d0 do 2105 ig=1,nalpha sum0=sum0+z(ig,t,i)*ralpha0(ig,j,k) 2105 continue gamz(j,k,t,i)=sum0 2104 continue c write(*,*)(gamz(j,ik,t,i),ik=1,K1-1) 2103 continue endif 2102 continue 2100 continue do 2106 i=1,N if(nr(i).ge.2) then do 2107 t=2,nr(i) do 2107 jl=1,K1 temp0=0.d0 do 2108 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,il,t,i)*yv(jl)+ & gamz(2,il,t,i)*yv(jl)**2.d0) 2108 continue do 2109 ik=1,K1 if(ik.lt.K1) then h0(jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,ik,t,i)*yv(jl)+ & gamz(2,ik,t,i)*yv(jl)**2.d0)/(1.d0+temp0) else h0(jl,K1,t,i)=1.d0/(1.d0+temp0) endif 2109 continue 2107 continue endif 2106 continue do 2150 i=1,N if(nr(i).ge.2) then do 2147 t=2,nr(i) do 2110 k=1,K1-1 do 2113 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2115 j=1,K1 if(kg.eq.k) dht0=h0(j,k,t,i)*(1.d0-h0(j,k,t,i)) if(kg.ne.k) dht0=-h0(j,k,t,i)*h0(j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0*pi(j,t-1,i) 2115 continue 2113 continue 2110 continue do 2120 k=1,K1-1 do 2120 ig=1,nx-1 psum(ig,k)=0.d0 2120 continue do 2123 k=1,K1-1 do 2125 j=1,K1 do 2130 ig=1,nx-1 if(j.eq.1) psum(ig,k)=psum(ig,k)+ & h0(j,k,t,i)*dpmdbet(ig,j,t-1,i) if((j.gt.1).and.(j.lt.K1)) psum(ig,k)=psum(ig,k)+ & h0(j,k,t,i)*(dpmdbet(ig,j,t-1,i)-dpmdbet(ig,j-1,t-1,i)) if(j.eq.K1) psum(ig,k)=psum(ig,k)+ & h0(j,k,t,i)*(-dpmdbet(ig,j-1,t-1,i)) 2130 continue 2125 continue do 2133 ig=1,nx-1 if((k.gt.1).and.(k.lt.K1)) psum(ig,k)=dpmdbet(ig,k,t,i)- & dpmdbet(ig,k-1,t,i)-psum(ig,k) if(k.eq.1) psum(ig,k)=dpmdbet(ig,k,t,i)-psum(ig,k) if(k.eq.K1) psum(ig,k)=-dpmdbet(ig,k-1,t,i)-psum(ig,k) 2133 continue 2123 continue do 2135 ik=1,(nx-1)*(K1-1) do 2135 ig=1,(nx-1)*(K1-1) dhtsum1(ik,ig)=0.d0 2135 continue do 2137 k=1,K1-1 do 2137 kg=1,K1-1 do 2137 ig=1,nx-1 dhtsum1((k-1)*(nx-1)+ig,(kg-1)*(nx-1)+ig)=dhtsum(k,kg) 2137 continue do 2140 ik=1,K1-1 do 2140 l=1,nx-1 psum1((ik-1)*(nx-1)+l)=psum(l,ik) 2140 continue call DLINRG((K1-1)*(nx-1),dhtsum1,(K1-1)*(nx-1), & Ainv,(K1-1)*(nx-1)) call DMURRV((K1-1)*(nx-1),(K1-1)*(nx-1),Ainv, & (K1-1)*(nx-1),(K1-1)*(nx-1),psum1,1,(K1-1)*(nx-1),res) do 2143 k=1,K1-1 do 2145 ig=1,nx-1 dtrdbet(ig,k,t,i)=res((k-1)*(nx-1)+ig) 2145 continue 2143 continue 2147 continue endif 2150 continue return end c Calculate d tri / d beta_0 subroutine dtrbet0(z,ralpha0,triangle,yv,pi,dpmdbet0,N,IT, & nr,K1,nalpha,ialph,dtrdbet0) integer t,nr(N) double precision z(nalpha,IT,N), & gamz(ialph,K1-1,IT,N),ralpha0(nalpha,ialph,K1-1), & triangle(K1-1,IT,N),yv(K1),h0(K1,K1,IT,N), & pi(K1,IT,N), & dpmdbet0(K1-1,K1-1,IT,N),dhtsum(K1-1,K1-1), & psum(K1-1),Ainv(K1-1,K1-1),resbet0(K1-1), & dtrdbet0(K1-1,K1-1,IT,N),sum0,temp0,dht0 do 2200 j=1,ialph do 2202 i=1,N if(nr(i).ge.2) then do 2203 t=2,nr(i) do 2204 k=1,K1-1 sum0=0.d0 do 2205 ig=1,nalpha sum0=sum0+z(ig,t,i)*ralpha0(ig,j,k) 2205 continue gamz(j,k,t,i)=sum0 2204 continue 2203 continue endif 2202 continue 2200 continue do 2206 i=1,N if(nr(i).ge.2) then do 2207 t=2,nr(i) do 2207 jl=1,K1 temp0=0.d0 do 2208 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,il,t,i)*yv(jl)+ & gamz(2,il,t,i)*yv(jl)**2.d0) 2208 continue do 2209 ik=1,K1 if(ik.lt.K1) then h0(jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,ik,t,i)*yv(jl)+ & gamz(2,ik,t,i)*yv(jl)**2.d0)/(1.d0+temp0) else h0(jl,K1,t,i)=1.d0/(1.d0+temp0) endif 2209 continue 2207 continue endif 2206 continue do 2210 i=1,N if(nr(i).ge.2) then do 2213 t=2,nr(i) do 2215 kg=1,K1-1 do 2217 k=1,K1-1 dhtsum(k,kg)=0.d0 do 2220 j=1,K1 if(kg.eq.k) dht0=h0(j,k,t,i)*(1.0-h0(j,k,t,i)) if(kg.ne.k) dht0=-h0(j,k,t,i)*h0(j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0*pi(j,t-1,i) 2220 continue 2217 continue 2215 continue c d tri / d beta_0 do 2225 ia=1,K1-1 do 2227 k=1,K1-1 psum(k)=0.d0 do 2230 j=1,K1 if(j.eq.1) psum(k)=psum(k)+h0(j,k,t,i)*dpmdbet0(ia,j,t-1,i) if((j.gt.1).and.(j.lt.K1)) psum(k)=psum(k)+ & h0(j,k,t,i)*(dpmdbet0(ia,j,t-1,i)-dpmdbet0(ia,j-1,t-1,i)) if(j.eq.K1) psum(k)=psum(k)+ & h0(j,k,t,i)*(-dpmdbet0(ia,j-1,t-1,i)) 2230 continue if((k.gt.1).and.(k.lt.K1)) psum(k)=dpmdbet0(ia,k,t,i)- & dpmdbet0(ia,k-1,t,i)-psum(k) if(k.eq.1) psum(k)=dpmdbet0(ia,k,t,i)-psum(k) if(k.eq.K1) psum(k)=-dpmdbet0(ia,k-1,t,i)-psum(k) 2227 continue call DLINRG(K1-1,dhtsum,K1-1,Ainv,K1-1) call DMURRV(K1-1,K1-1,Ainv,K1-1,K1-1,psum,1,K1-1,resbet0) do 2235 ik=1,K1-1 dtrdbet0(ia,ik,t,i)=resbet0(ik) 2235 continue 2225 continue 2213 continue endif 2210 continue return end c Calculate d tri / d alpha_1^(a) subroutine dtralp1(z,ralpha0,triangle,yv,pi,dhdgam1, & N,IT,nr,K1,nalpha,ialph,dtrdalp1) integer t,nr(N) double precision z(nalpha,IT,N),yv(K1), & ralpha0(nalpha,ialph,K1-1), & h0(K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,K1-1,IT,N), & triangle(K1-1,IT,N),dhdgam1(K1-1,K1,K1,IT,N), & dhtsum(K1-1,K1-1),dtrdalp1(nalpha,K1-1,K1-1,IT,N), & psum(nalpha,K1-1), & Ainv((K1-1)*nalpha,(K1-1)*nalpha),res((K1-1)*nalpha), & dhtsum1(nalpha*(K1-1),nalpha*(K1-1)), & psum1((K1-1)*nalpha),sum0,temp0,dht0 do 2300 j=1,ialph do 2302 i=1,N if(nr(i).ge.2) then do 2303 t=2,nr(i) do 2304 k=1,K1-1 sum0=0.d0 do 2305 ig=1,nalpha sum0=sum0+z(ig,t,i)*ralpha0(ig,j,k) 2305 continue gamz(j,k,t,i)=sum0 2304 continue 2303 continue endif 2302 continue 2300 continue do 2306 i=1,N if(nr(i).ge.2) then do 2307 t=2,nr(i) do 2307 jl=1,K1 temp0=0.d0 do 2308 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,il,t,i)*yv(jl)+ & gamz(2,il,t,i)*yv(jl)**2.d0) 2308 continue do 2309 ik=1,K1 if(ik.lt.K1) then h0(jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,ik,t,i)*yv(jl)+ & gamz(2,ik,t,i)*yv(jl)**2)/(1.0+temp0) else h0(jl,K1,t,i)=1.0/(1.0+temp0) endif 2309 continue 2307 continue endif 2306 continue do 2350 i=1,N if(nr(i).ge.2) then do 2347 t=2,nr(i) do 2310 k=1,K1-1 do 2313 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2315 j=1,K1 if(kg.eq.k) dht0=h0(j,k,t,i)*(1.0-h0(j,k,t,i)) if(kg.ne.k) dht0=-h0(j,k,t,i)*h0(j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0*pi(j,t-1,i) 2315 continue 2313 continue 2310 continue do 2317 ia=1,K1-1 do 2320 k=1,K1-1 do 2320 ig=1,nalpha psum(ig,k)=0.d0 2320 continue do 2323 k=1,K1-1 do 2330 in=1,nalpha do 2325 ig=1,K1 psum(in,k)=psum(in,k)- & dhdgam1(ia,ig,k,t,i)*z(in,t,i)*pi(ig,t-1,i) 2325 continue 2330 continue 2323 continue do 2335 ik=1,nalpha*(K1-1) do 2335 ig=1,nalpha*(K1-1) dhtsum1(ik,ig)=0.d0 2335 continue do 2337 k=1,K1-1 do 2337 kg=1,K1-1 do 2337 ig=1,nalpha dhtsum1((k-1)*nalpha+ig,(kg-1)*nalpha+ig)=dhtsum(k,kg) 2337 continue do 2340 ik=1,K1-1 do 2340 l=1,nalpha psum1((ik-1)*nalpha+l)=psum(l,ik) 2340 continue call DLINRG((K1-1)*nalpha,dhtsum1,(K1-1)*nalpha, & Ainv,(K1-1)*nalpha) call DMURRV((K1-1)*nalpha,(K1-1)*nalpha,Ainv, & (K1-1)*nalpha,(K1-1)*nalpha,psum1,1,(K1-1)*nalpha,res) do 2343 k=1,K1-1 do 2345 ig=1,nalpha dtrdalp1(ig,ia,k,t,i)=res((k-1)*nalpha+ig) 2345 continue 2343 continue 2317 continue 2347 continue endif 2350 continue return end c Calculate d tri / d alpha_2^(a) subroutine dtralp2(z,ralpha0,triangle,yv,pi,dhdgam2, & N,IT,nr,K1,nalpha,ialph,dtrdalp2) integer t,nr(N) double precision z(nalpha,IT,N),yv(K1), & ralpha0(nalpha,ialph,K1-1), & h0(K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,K1-1,IT,N), & triangle(K1-1,IT,N),dhdgam2(K1-1,K1,K1,IT,N), & dhtsum(K1-1,K1-1),dtrdalp2(nalpha,K1-1,K1-1,IT,N), & psum(nalpha,K1-1), & Ainv((K1-1)*nalpha,(K1-1)*nalpha),res((K1-1)*nalpha), & dhtsum1(nalpha*(K1-1),nalpha*(K1-1)), & psum1((K1-1)*nalpha),sum0,temp0,dht0 do 2400 j=1,ialph do 2402 i=1,N if(nr(i).ge.2) then do 2403 t=2,nr(i) do 2404 k=1,K1-1 sum0=0.d0 do 2405 ig=1,nalpha sum0=sum0+z(ig,t,i)*ralpha0(ig,j,k) 2405 continue gamz(j,k,t,i)=sum0 2404 continue 2403 continue endif 2402 continue 2400 continue do 2406 i=1,N if(nr(i).ge.2) then do 2407 t=2,nr(i) do 2407 jl=1,K1 temp0=0.d0 do 2408 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,il,t,i)*yv(jl)+ & gamz(2,il,t,i)*yv(jl)**2.0) 2408 continue do 2409 ik=1,K1 if(ik.lt.K1) then h0(jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,ik,t,i)*yv(jl)+ & gamz(2,ik,t,i)*yv(jl)**2)/(1.d0+temp0) else h0(jl,K1,t,i)=1.d0/(1.d0+temp0) endif 2409 continue 2407 continue endif 2406 continue do 2450 i=1,N if(nr(i).ge.2) then do 2447 t=2,nr(i) do 2410 k=1,K1-1 do 2413 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2415 j=1,K1 if(kg.eq.k) dht0=h0(j,k,t,i)*(1.d0-h0(j,k,t,i)) if(kg.ne.k) dht0=-h0(j,k,t,i)*h0(j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0*pi(j,t-1,i) 2415 continue 2413 continue 2410 continue do 2417 ib=1,K1-1 do 2420 k=1,K1-1 do 2420 ig=1,nalpha psum(ig,k)=0.d0 2420 continue do 2423 k=1,K1-1 do 2425 in=1,nalpha do 2430 ig=1,K1 psum(in,k)=psum(in,k)- & dhdgam2(ib,ig,k,t,i)*z(in,t,i)*pi(ig,t-1,i) 2430 continue 2425 continue 2423 continue do 2435 ik=1,nalpha*(K1-1) do 2435 ig=1,nalpha*(K1-1) dhtsum1(ik,ig)=0.d0 2435 continue do 2437 k=1,K1-1 do 2437 kg=1,K1-1 do 2437 ig=1,nalpha dhtsum1((k-1)*nalpha+ig,(kg-1)*nalpha+ig)=dhtsum(k,kg) 2437 continue do 2440 ik=1,K1-1 do 2440 l=1,nalpha psum1((ik-1)*nalpha+l)=psum(l,ik) 2440 continue call DLINRG((K1-1)*nalpha,dhtsum1,(K1-1)*nalpha, & Ainv,(K1-1)*nalpha) call DMURRV((K1-1)*nalpha,(K1-1)*nalpha,Ainv, & (K1-1)*nalpha,(K1-1)*nalpha,psum1,1,(K1-1)*nalpha,res) do 2443 k=1,K1-1 do 2445 ig=1,nalpha dtrdalp2(ig,ib,k,t,i)=res((k-1)*nalpha+ig) 2445 continue 2443 continue 2417 continue 2447 continue endif 2450 continue return end