c Generate IOMTM(2) for ordinal data c Fit IOMTM(2) c Fisher-scoring algorithm with linear and quadratic c Ignore orthogonality for Hessain matrix c i.e. Consider H with \beta_0, beta, alpha c Covariates : Visit and Arm c Real Ordinal data analysis $Debug use msimsl parameter (IT=8,K1=4,ialph=2,norder=2,M=1000,M1=5000, & N=500,nx=4,nalpha=2,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), & ba2((K1-1)*nalpha1*ialph*norder), & bb0(K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & 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.5d0,0.1d0, & -0.5d0 / data ba0 / -0.1d0,-0.3d0,0.5d0,-0.1d0,0.1d0,0.2d0, & -0.1d0,-0.1d0,0.3d0,-0.2d0,0.2d0,0.1d0, & -0.1d0,-0.2d0,0.5d0,-0.1d0,0.3d0,0.1d0, & -0.1d0,-0.4d0,0.3d0,-0.2d0,0.2d0,0.2d0 / data ba2 / -1.8d0,-1.1d0,-0.47d0, & 0.07d0,-0.4d0,-0.47d0, & -1.8d0,-1.1d0,-0.47d0, & 0.07d0,-0.4d0,-0.47d0 / 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 x(2,t,i)=(dble(t)-1.d0)/10.d0 c Arm transformation (IFL=control group, FOLFOX, IROX=treatment group) if(i.lt.(N/3))then c group=1 x(3,t,i)=0.d0 x(4,t,i)=0.d0 else if(i.lt.(N/3*2))then c group=2 x(3,t,i)=1.d0 x(4,t,i)=0.d0 else c group=3 x(3,t,i)=0.d0 x(4,t,i)=1.d0 endif 23 continue 20 continue open(3,file='Fit-MTM2_Sim-MTM2_MAR.out') do 70 iter=1,M 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 write(*,*)'iter=',iter do 25 i=1,N nr(i)=IT do 25 t=1,nr(i) z(1,t,i)=1.d0 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 k=1,K1-1 do 27 j=1,ialph 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)=b1(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 caltri0(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.2d0*y1(t-1,i)+0.1d0*y1(t-2,i))/ & (1.d0+dexp(-2.5d0+0.2d0*y1(t-1,i)+0.1d0*y1(t-2,i))) pdrop(2)=1.d0/ & (1.d0+dexp(-2.5d0+0.2d0*y1(t-1,i)+0.1d0*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) enddo write(*,*)inr c read(*,*)iii call findmle(x,y,yv,z,r,id,nr,b1,ba2,bb0, & N,IT,K1,nx,nalpha1,ialph,norder,ifail) c call findmle(x,y,yv,z,r,id,nr,b1,ba1,bb0, c & N,IT,K1,nx,nalpha,ialph,norder,ifail) if(ifail.eq.0) then write(*,75)(bb0(ig),ig=1,K1-1+nx-1) write(3,75)(bb0(ig),ig=1,K1-1+nx-1) 75 format(100f12.4) endif write(*,*)'After findmle' write(*,78)(ba1(ig),ig=1,(K1-1)*nalpha*ialph*norder) 78 format(100f12.4) 70 continue close(3) stop end c Subroutine for calculating \triangle_{itk} c for findcand subroutine caltri0(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----------------------------------------------------------c c Find Mle for beta_0, beta, alpha c c----------------------------------------------------------c subroutine findmle(x,y,yv,z,r,id,nr,b1,ba1,bb1, & N,IT,K1,nx,nalpha,ialph,norder,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), & b1(K1-1+nx-1),ba1((K1-1)*nalpha*ialph*norder), & bb0(K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & bb1(K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & ralpha0(nalpha,ialph,K1-1,norder), & gamz(ialph,norder,K1-1,IT,N), & triangle(K1-1,IT,N),pm(K1,IT,N),pi(K1,IT,N), & beta1(nx,K1), & 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,ialph,K1-1,K1-1,IT,N), & dtrdalp2(nalpha,ialph,K1-1,K1-1,IT,N), & dhdgam1(ialph,K1-1,K1,K1,K1,IT,N), & dhdgam2(ialph,K1-1,K1,K1,K1,IT,N), & sumalp1(nalpha,ialph,K1-1),sumalp2(nalpha,ialph,K1-1), & dQ1(K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & H1(K1-1+nx-1+(K1-1)*nalpha*ialph*norder, & K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & H1inv(K1-1+nx-1+(K1-1)*nalpha*ialph*norder, & K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & tempb1(K1-1+nx-1+(K1-1)*nalpha*ialph*norder), & pj(K1,K1,IT,N), & pc(K1,IT,N),h(K1,K1,K1,IT,N),dh(K1-1,K1,K1,K1,IT,N), & dl(K1-1+nx-1),ddb1(nx-1,nx-1), & ddb0(K1-1,K1-1),ddb10(nx-1,K1-1), & sumexph(K1,K1),sumbet(nx-1),sumbet0(K1-1), & sumexp,tempdh,temp0,temp1,temp2,aloglik,aloglik0 do 480 ig=1,K1-1+nx-1 bb0(ig)=b1(ig) 480 continue do 481 ig=1,(K1-1)*nalpha*ialph*norder bb0(K1-1+nx-1+ig)=ba1(ig) 481 continue do 482 l=1,norder do 482 k=1,K1-1 do 482 j=1,ialph do 482 ig=1,nalpha ralpha0(ig,j,k,l)= & ba1((K1-1)*ialph*nalpha*(l-1)+ialph*nalpha*(k-1)+ & (j-1)*nalpha+ig) 482 continue c write(*,*)'beta1=' 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 c write(*,*)(beta1(ig,k),ig=1,nx) 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,pj,ralpha0,gamz,N,nr,IT,K1,ialph,nalpha, & norder,triangle,2) do 509 i=1,N do 509 t=1,IT do 509 k=1,K1 do 509 j1=1,K1 do 509 j2=1,K1 h(j2,j1,k,t,i)=0.d0 509 continue c Calculate P_{itk}^c aloglik0=aloglik aloglik=0.d0 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) if(t.eq.2) then sumexp=0.d0 do 532 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) 532 continue do 533 k=1,K1-1 pc(k,t,i)=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) 533 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)+gamz(1,1,ik1,t,i)*y(t-1,i)+ & gamz(2,1,ik1,t,i)*y(t-1,i)**2.d0) 534 continue aloglik=aloglik+dlog(pc(K1,t,i)) do 535 j=1,K1 sumexph(1,j)=0.d0 do 536 l=1,K1-1 sumexph(1,j)=sumexph(1,j)+dexp(triangle(l,t,i)+ & gamz(1,1,l,t,i)*yv(j)+gamz(2,1,l,t,i)*yv(j)**2.d0) 536 continue do 537 k=1,K1-1 h(1,j,k,t,i)=dexp(triangle(k,t,i)+ & gamz(1,1,k,t,i)*yv(j)+gamz(2,1,k,t,i)*yv(j)**2.d0)/ & (1.d0+sumexph(1,j)) 537 continue h(1,j,K1,t,i)=1.d0/(1.d0+sumexph(1,j)) 535 continue else sumexp=0.d0 do 538 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) 538 continue do 539 k=1,K1-1 pc(k,t,i)=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) 539 continue pc(K1,t,i)=1.d0/(1.d0+sumexp) do 540 ik1=1,K1-1 aloglik=aloglik+dble(id(ik1,t,i))* & (triangle(ik1,t,i)+gamz(1,1,ik1,t,i)*y(t-1,i)+ & gamz(2,1,ik1,t,i)*y(t-1,i)**2.d0+ & gamz(1,2,ik1,t,i)*y(t-2,i)+ & gamz(2,2,ik1,t,i)*y(t-2,i)**2.d0) 540 continue aloglik=aloglik+dlog(pc(K1,t,i)) do 541 j1=1,K1 do 541 j2=1,K1 sumexph(j2,j1)=0.d0 do 542 l=1,K1-1 sumexph(j2,j1)=sumexph(j2,j1)+dexp(triangle(l,t,i)+ & 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) 542 continue do 543 k=1,K1-1 h(j2,j1,k,t,i)=dexp(triangle(k,t,i)+ & gamz(1,1,k,t,i)*yv(j1)+gamz(2,1,k,t,i)*yv(j1)**2.d0+ & gamz(1,2,k,t,i)*yv(j2)+gamz(2,2,k,t,i)*yv(j2)**2.d0)/ & (1.d0+sumexph(j2,j1)) 543 continue h(j2,j1,K1,t,i)=1.d0/(1.d0+sumexph(j2,j1)) 541 continue endif 530 continue endif 527 continue write(*,*)'loglikelihood =',aloglik if(aloglik.lt.aloglik0) then nolog=nolog+1 if(nolog.gt.20) then ifail=1 goto 860 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.d0-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_{ikgl}^(t) / d gamma_{it1j}^(a) and d h_{ikgl}^(t) / d gamma_{it2j}^(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 ig1=1,K1 do 549 ig2=1,K1 do 549 ia=1,K1-1 do 549 ib=1,ialph if(t.eq.2) then dhdgam1(ib,ia,ig2,ig1,k,t,i)=0.d0 dhdgam2(ib,ia,ig2,ig1,k,t,i)=0.d0 if(k.eq.ia)then dhdgam1(ib,ia,1,ig1,k,t,i)=h(1,ig1,k,t,i)* & (1.d0-h(1,ig1,k,t,i))*yv(ig1)**dble(ib) else dhdgam1(ib,ia,1,ig1,k,t,i)=-h(1,ig1,k,t,i)* & h(1,ig1,ia,t,i)*yv(ig1)**dble(ib) endif else if(k.eq.ia)then dhdgam1(ib,ia,ig2,ig1,k,t,i)=h(ig2,ig1,k,t,i)* & (1.d0-h(ig2,ig1,k,t,i))*yv(ig1)**dble(ib) dhdgam2(ib,ia,ig2,ig1,k,t,i)=h(ig2,ig1,k,t,i)* & (1.d0-h(ig2,ig1,k,t,i))*yv(ig2)**dble(ib) else dhdgam1(ib,ia,ig2,ig1,k,t,i)=-h(ig2,ig1,k,t,i)* & h(ig2,ig1,ia,t,i)*yv(ig1)**dble(ib) dhdgam2(ib,ia,ig2,ig1,k,t,i)=-h(ig2,ig1,k,t,i)* & h(ig2,ig1,ia,t,i)*yv(ig2)**dble(ib) endif endif 549 continue endif 548 continue c Calculate d tri / d beta call dtrbet(gamz,triangle,yv,pi,pj,dpmdbet, & N,IT,nr,K1,norder,nx,ialph,dtrdbet) c Calculate d tri / d beta_0 call dtrbet0(gamz,triangle,yv,pi,pj,dpmdbet0,N,IT, & nr,K1,norder,ialph,dtrdbet0) c Calculate d tri / d alpha_1^(a) call dtralp1(z,triangle,yv,pi,pj,dhdgam1,gamz,N,IT, & nr,K1,nalpha,ialph,norder,dtrdalp1) c Calculate d tri / d alpha_2^(b) call dtralp2(z,triangle,yv,pi,pj,dhdgam2,gamz,N,IT, & nr,K1,nalpha,ialph,norder,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 ib=1,ialph do 552 j=1,nalpha sumalp1(j,ib,ia)=0.d0 sumalp2(j,ib,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 556 k=1,K1-1 sumbet0(j)=sumbet0(j)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdbet0(j,k,t,i) 556 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 ib=1,ialph 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,ib,ia)=sumalp1(ig,ib,ia)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdalp1(ig,ib,ia,k,t,i) if(ia.eq.k) then sumalp1(ig,ib,ia)=sumalp1(ig,ib,ia)+ & (dble(id(ia,t,i))-pc(ia,t,i))*y(t-1,i)**dble(ib)*z(ig,t,i) endif if(t.gt.2) then sumalp2(ig,ib,ia)=sumalp2(ig,ib,ia)+ & (dble(id(k,t,i))-pc(k,t,i))*dtrdalp2(ig,ib,ia,k,t,i) if(ia.eq.k) then sumalp2(ig,ib,ia)=sumalp2(ig,ib,ia)+ & (dble(id(ia,t,i))-pc(ia,t,i))*y(t-2,i)**dble(ib)*z(ig,t,i) endif endif 567 continue 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 ib=1,ialph do 570 ig=1,nalpha dQ1(K1-1+nx-1+(ia-1)*ialph*nalpha+(ib-1)*nalpha+ig)= & sumalp1(ig,ib,ia) dQ1(K1-1+nx-1+(K1-1)*ialph*nalpha+(ia-1)*ialph*nalpha+ & (ib-1)*nalpha+ig)=sumalp2(ig,ib,ia) 570 continue c Calculate Hessian matrix do 571 ia1=1,K1-1+nx-1+norder*(K1-1)*ialph*nalpha do 571 ia2=1,K1-1+nx-1+norder*(K1-1)*ialph*nalpha H1(ia1,ia2)=0.d0 571 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 j1=1,K1 do 585 j2=1,K1 if(t.eq.2) then dh(ig,j2,j1,k,t,i)=0.d0 if(ig.eq.k) then dh(ig,1,j1,k,t,i)=h(1,j1,k,t,i)*(1.d0-h(1,j1,k,t,i)) else if(ig.ne.k) then dh(ig,1,j1,k,t,i)=-h(1,j1,k,t,i)*h(1,j1,ig,t,i) endif else if(ig.eq.k) then dh(ig,j2,j1,k,t,i)=h(j2,j1,k,t,i)*(1.d0-h(j2,j1,k,t,i)) else if(ig.ne.k) then dh(ig,j2,j1,k,t,i)=-h(j2,j1,k,t,i)*h(j2,j1,ig,t,i) endif 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 if(t.eq.2) then do 595 ik=1,K1 tempdh=tempdh+dh(ig,1,ik,k,t,i)*pi(ik,t-1,i) 595 continue else do 596 ik1=1,K1 do 596 ik2=1,K1 tempdh=tempdh+dh(ig,ik2,ik1,k,t,i)*pj(ik2,ik1,t-1,i) 596 continue endif c write(*,*)'tempdh=',tempdh H1(i1,i2)=H1(i1,i2)+ & tempdh*dtrdbet0(i1,ig,t,i)*dtrdbet0(i2,k,t,i) 594 continue 593 continue 591 continue endif 590 continue 587 continue do 598 i1=1,K1-1 do 598 i2=1,K1-1 H1(i1,i2)=H1(i1,i2)+ddb0(i1,i2) 598 continue c write(*,*)'sumbet0' c write(*,*)sumbet0 c write(*,*)'Hessian matrix for beta0' c do ig1=1,K1-1 c write(*,6000)(H1(ig1,ig2),ig2=1,K1-1) c6000 format(10f20.5) c enddo c \beta \beta do 600 i1=1,nx-1 do 600 i2=1,nx-1 do 602 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 if(t.eq.2) then do 608 ik=1,K1 tempdh=tempdh+dh(ig,1,ik,k,t,i)*pi(ik,t-1,i) 608 continue else do 609 ik1=1,K1 do 609 ik2=1,K1 tempdh=tempdh+dh(ig,ik2,ik1,k,t,i)*pj(ik2,ik1,t-1,i) 609 continue endif H1(K1-1+i1,K1-1+i2)=H1(K1-1+i1,K1-1+i2)+ & tempdh*dtrdbet(i1,ig,t,i)*dtrdbet(i2,k,t,i) 607 continue 605 continue 603 continue endif 602 continue 600 continue do 599 i1=1,nx-1 do 599 i2=1,nx-1 H1(K1-1+i1,K1-1+i2)=H1(K1-1+i1,K1-1+i2)+ddb1(i1,i2) 599 continue c write(*,*)'beta' c do ig1=1,nx-1 c write(*,6001)(H1(K1-1+ig1,K1-1+ig2),ig2=1,nx-1) c6001 format(10f15.5) c enddo 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 if(t.eq.2) then do 620 ik=1,K1 tempdh=tempdh+dh(ig,1,ik,k,t,i)*pi(ik,t-1,i) 620 continue else do 621 ik1=1,K1 do 621 ik2=1,K1 tempdh=tempdh+dh(ig,ik2,ik1,k,t,i)*pj(ik2,ik1,t-1,i) 621 continue endif H1(K1-1+i1,i2)=H1(K1-1+i1,i2)+ & tempdh*dtrdbet(i1,ig,t,i)*dtrdbet0(i2,k,t,i) 617 continue 615 continue 613 continue endif 611 continue 610 continue do 622 i1=1,nx-1 do 622 i2=1,K1-1 H1(K1-1+i1,i2)=H1(K1-1+i1,i2)+ddb10(i1,i2) 622 continue c \alpha_1b^(a) \alph_1b'^(a') do 623 ia1=1,K1-1 do 623 ia2=1,K1-1 do 623 ib1=1,ialph do 623 ib2=1,ialph 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 if(t.eq.2) then do 632 ij=1,K1 tempdh=tempdh+dh(ig,1,ij,k,t,i)*pi(ij,t-1,i) 632 continue else do 633 ij1=1,K1 do 633 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 633 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & tempdh*dtrdalp1(ig1,ib1,ia1,ig,t,i)* & dtrdalp1(ig2,ib2,ia2,k,t,i) 630 continue temp0=0.d0 if(t.eq.2) then do 635 ij=1,K1 temp0=temp0+dhdgam1(ib1,ia1,1,ij,k,t,i)*pi(ij,t-1,i) 635 continue else do 636 ij1=1,K1 do 636 ij2=1,K1 temp0=temp0+dhdgam1(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 636 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp0*z(ig1,t,i)*dtrdalp1(ig2,ib2,ia2,k,t,i) 628 continue do 637 ig=1,K1-1 temp1=0.d0 if(t.eq.2) then do 640 ij=1,K1 temp1=temp1+dh(ig,1,ij,ia2,t,i)* & yv(ij)**dble(ib2)*pi(ij,t-1,i) 640 continue else do 641 ij1=1,K1 do 641 ij2=1,K1 temp1=temp1+dh(ig,ij2,ij1,ia2,t,i)* & yv(ij1)**dble(ib2)*pj(ij2,ij1,t-1,i) 641 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp1*dtrdalp1(ig1,ib1,ia1,ig,t,i)*z(ig2,t,i) 637 continue temp2=0.d0 if(t.eq.2) then do 643 ij=1,K1 temp2=temp2+dhdgam1(ib1,ia1,1,ij,ia2,t,i)* & yv(ij)**dble(ib2)*pi(ij,t-1,i) 643 continue else do 644 ij1=1,K1 do 644 ij2=1,K1 temp2=temp2+dhdgam1(ib1,ia1,ij2,ij1,ia2,t,i)* & yv(ij1)**dble(ib2)*pj(ij2,ij1,t-1,i) 644 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp2*z(ig1,t,i)*z(ig2,t,i) 626 continue endif 625 continue 623 continue c \alpha_2 \alph_2 do 650 ia1=1,K1-1 do 650 ia2=1,K1-1 do 650 ib1=1,ialph do 650 ib2=1,ialph do 650 ig1=1,nalpha do 650 ig2=1,nalpha do 652 i=1,N if(nr(i).ge.3) then do 654 t=3,nr(i) do 658 k=1,K1-1 do 660 ig=1,K1-1 tempdh=0.d0 do 663 ij1=1,K1 do 663 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 663 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & tempdh*dtrdalp2(ig1,ib1,ia1,ig,t,i)* & dtrdalp2(ig2,ib2,ia2,k,t,i) 660 continue temp0=0.d0 do 666 ij1=1,K1 do 666 ij2=1,K1 temp0=temp0+dhdgam2(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 666 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp0*z(ig1,t,i)*dtrdalp2(ig2,ib2,ia2,k,t,i) 658 continue do 667 ig=1,K1-1 temp1=0.d0 do 671 ij1=1,K1 do 671 ij2=1,K1 temp1=temp1+dh(ig,ij2,ij1,ia2,t,i)* & yv(ij2)**dble(ib2)*pj(ij2,ij1,t-1,i) 671 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp1*dtrdalp2(ig1,ib1,ia1,ig,t,i)*z(ig2,t,i) 667 continue temp2=0.d0 do 674 ij1=1,K1 do 674 ij2=1,K1 temp2=temp2+dhdgam2(ib1,ia1,ij2,ij1,ia2,t,i)* & yv(ij2)**dble(ib2)*pj(ij2,ij1,t-1,i) 674 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp2*z(ig1,t,i)*z(ig2,t,i) 654 continue endif 652 continue 650 continue c \alpha_2 \alph_1 do 680 ia1=1,K1-1 do 680 ia2=1,K1-1 do 680 ib1=1,ialph do 680 ib2=1,ialph do 680 ig1=1,nalpha do 680 ig2=1,nalpha do 682 i=1,N if(nr(i).ge.3) then do 684 t=3,nr(i) do 688 k=1,K1-1 do 690 ig=1,K1-1 tempdh=0.d0 do 693 ij1=1,K1 do 693 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 693 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & tempdh*dtrdalp2(ig1,ib1,ia1,ig,t,i)* & dtrdalp1(ig2,ib2,ia2,k,t,i) 690 continue temp0=0.d0 do 696 ij1=1,K1 do 696 ij2=1,K1 temp0=temp0+dhdgam2(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 696 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp0*z(ig1,t,i)*dtrdalp1(ig2,ib2,ia2,k,t,i) 688 continue do 697 ig=1,K1-1 temp1=0.d0 do 701 ij1=1,K1 do 701 ij2=1,K1 temp1=temp1+dh(ig,ij2,ij1,ia2,t,i)* & yv(ij1)**dble(ib2)*pj(ij2,ij1,t-1,i) 701 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp1*dtrdalp2(ig1,ib1,ia1,ig,t,i)*z(ig2,t,i) 697 continue temp2=0.d0 do 704 ij1=1,K1 do 704 ij2=1,K1 temp2=temp2+dhdgam2(ib1,ia1,ij2,ij1,ia2,t,i)* & yv(ij1)**dble(ib2)*pj(ij2,ij1,t-1,i) 704 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+nx-1+(ia2-1)*ialph*nalpha+(ib2-1)*nalpha+ig2)+ & temp2*z(ig1,t,i)*z(ig2,t,i) 684 continue endif 682 continue 680 continue c alpha_1 beta do 720 ia1=1,K1-1 do 720 ib1=1,ialph do 720 ig1=1,nalpha do 720 ig2=1,nx-1 do 722 i=1,N if(nr(i).ge.2) then do 724 t=2,nr(i) do 726 k=1,K1-1 do 728 ig=1,K1-1 tempdh=0.d0 if(t.eq.2) then do 730 ij=1,K1 tempdh=tempdh+dh(ig,1,ij,k,t,i)*pi(ij,t-1,i) 730 continue else do 732 ij1=1,K1 do 732 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 732 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+ig2)+ & tempdh*dtrdalp1(ig1,ib1,ia1,ig,t,i)* & dtrdbet(ig2,k,t,i) 728 continue temp0=0.d0 if(t.eq.2) then do 735 ij=1,K1 temp0=temp0+dhdgam1(ib1,ia1,1,ij,k,t,i)*pi(ij,t-1,i) 735 continue else do 736 ij1=1,K1 do 736 ij2=1,K1 temp0=temp0+dhdgam1(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 736 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & K1-1+ig2)+ & temp0*z(ig1,t,i)*dtrdbet(ig2,k,t,i) 726 continue 724 continue endif 722 continue 720 continue c alpha_2 beta do 740 ia1=1,K1-1 do 740 ib1=1,ialph do 740 ig1=1,nalpha do 740 ig2=1,nx-1 do 742 i=1,N if(nr(i).ge.3) then do 744 t=3,nr(i) do 746 k=1,K1-1 do 748 ig=1,K1-1 tempdh=0.d0 do 752 ij1=1,K1 do 752 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 752 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,K1-1+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,K1-1+ig2)+ & tempdh*dtrdalp2(ig1,ib1,ia1,ig,t,i)* & dtrdbet(ig2,k,t,i) 748 continue temp0=0.d0 do 756 ij1=1,K1 do 756 ij2=1,K1 temp0=temp0+dhdgam2(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 756 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,K1-1+ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,K1-1+ig2)+ & temp0*z(ig1,t,i)*dtrdbet(ig2,k,t,i) 746 continue 744 continue endif 742 continue 740 continue c alpha_1 beta_0j do 760 ia1=1,K1-1 do 760 ib1=1,ialph do 760 ig1=1,nalpha do 760 ig2=1,K1-1 do 762 i=1,N if(nr(i).ge.2) then do 764 t=2,nr(i) do 766 k=1,K1-1 do 768 ig=1,K1-1 tempdh=0.d0 if(t.eq.2) then do 770 ij=1,K1 tempdh=tempdh+dh(ig,1,ij,k,t,i)*pi(ij,t-1,i) 770 continue else do 772 ij1=1,K1 do 772 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 772 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1, & ig2)+ & tempdh*dtrdalp1(ig1,ib1,ia1,ig,t,i)* & dtrdbet0(ig2,k,t,i) 768 continue temp0=0.d0 if(t.eq.2) then do 775 ij=1,K1 temp0=temp0+dhdgam1(ib1,ia1,1,ij,k,t,i)*pi(ij,t-1,i) 775 continue else do 776 ij1=1,K1 do 776 ij2=1,K1 temp0=temp0+dhdgam1(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 776 continue endif H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,ig2)= & H1(K1-1+nx-1+(ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,ig2)+ & temp0*z(ig1,t,i)*dtrdbet0(ig2,k,t,i) 766 continue 764 continue endif 762 continue 760 continue c alpha_2 beta_0j do 780 ia1=1,K1-1 do 780 ib1=1,ialph do 780 ig1=1,nalpha do 780 ig2=1,K1-1 do 782 i=1,N if(nr(i).ge.3) then do 784 t=3,nr(i) do 786 k=1,K1-1 do 788 ig=1,K1-1 tempdh=0.d0 do 792 ij1=1,K1 do 792 ij2=1,K1 tempdh=tempdh+dh(ig,ij2,ij1,k,t,i)*pj(ij2,ij1,t-1,i) 792 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,ig2)= * H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,ig2)+ & tempdh*dtrdalp2(ig1,ib1,ia1,ig,t,i)* & dtrdbet0(ig2,k,t,i) 788 continue temp0=0.d0 do 796 ij1=1,K1 do 796 ij2=1,K1 temp0=temp0+dhdgam2(ib1,ia1,ij2,ij1,k,t,i)* & pj(ij2,ij1,t-1,i) 796 continue H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,ig2)= & H1(K1-1+nx-1+(K1-1)*ialph*nalpha+ & (ia1-1)*ialph*nalpha+(ib1-1)*nalpha+ig1,ig2)+ & temp0*z(ig1,t,i)*dtrdbet0(ig2,k,t,i) 786 continue 784 continue endif 782 continue 780 continue do 800 ig1=1,K1-1+nx-1+(K1-1)*ialph*nalpha*norder do 800 ig2=1,ig1 H1(ig2,ig1)=H1(ig1,ig2) 800 continue call DLINRG(K1-1+nx-1+(K1-1)*ialph*nalpha*norder,H1, & K1-1+nx-1+(K1-1)*ialph*nalpha*norder,H1inv, & K1-1+nx-1+(K1-1)*ialph*nalpha*norder) call DMURRV(K1-1+nx-1+(K1-1)*ialph*nalpha*norder, & K1-1+nx-1+(K1-1)*ialph*nalpha*norder,H1inv, & K1-1+nx-1+(K1-1)*ialph*nalpha*norder, & K1-1+nx-1+(K1-1)*ialph*nalpha*norder,dQ1,1, & K1-1+nx-1+(K1-1)*ialph*nalpha*norder,tempb1) do 833 ig=1,K1-1+nx-1+(K1-1)*ialph*nalpha*norder bb1(ig)=bb0(ig)+tempb1(ig) 833 continue c write(*,*)'Derivatives of log likelihoo for marginal' c write(*,*)(dQ1(ig),ig=1,K1-1+nx-1) c write(*,*)'Derivatives of log likelihoo for conditional' c write(*,*) c & (dQ1(ig),ig=K1-1+nx-1+1,K1-1+nx-1+(K1-1)*nalpha*ialph*norder) c write(*,*)'Estimate beta' c write(*,*)(bb1(ig),ig=1,K1-1+nx-1) c write(*,*)'Estimate alpha' c write(*,*)(bb1(K1-1+nx-1+ig),ig=1,(K1-1)*nalpha*ialph*norder) update=0.d0 do 835 ig=1,K1-1+nx-1+(K1-1)*nalpha*ialph*norder update=update+abs(bb1(ig)-bb0(ig)) 835 continue c write(*,*)'update=',update c write(*,*)'------------------------------' if(update.le.0.0001d0) then ifail=0 goto 850 endif do 837 ig=1,K1-1+nx-1 b1(ig)=bb1(ig) 837 continue do 848 ig=1,(K1-1)*nalpha*ialph*norder ba1(ig)=bb1(K1-1+nx-1+ig) 848 continue do 849 ig=1,K1-1+nx-1+(K1-1)*nalpha*ialph*norder bb0(ig)=bb1(ig) 849 continue do 840 l=1,norder do 840 k=1,K1-1 do 840 j=1,ialph do 840 ig=1,nalpha ralpha0(ig,j,k,l)= & ba1((K1-1)*ialph*nalpha*(l-1)+ialph*nalpha*(k-1)+ & (j-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' open(2,file='hhh.out') do ig=1,K1-1+nx-1+(K1-1)*nalpha*ialph*norder write(2,851)(H1(ig,ik), & ik=1,K1-1+nx-1+(K1-1)*nalpha*ialph*norder) 851 format(100f10.3) enddo close(2) 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 dl(ig)=dlike0(ig) 2090 continue do 2093 ig=K1,nx-1+K1-1 dl(ig)=dlike1(ig-(K1-1)) 2093 continue return end c Calculate d tri / d beta subroutine dtrbet(gamz,triangle,yv,pi,pj,dpmdbet, & N,IT,nr,K1,norder,nx,ialph,dtrdbet) integer t,nr(N) double precision pj(K1,K1,IT,N), & h(K1,K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,norder,K1-1,IT,N), & triangle(K1-1,IT,N),yv(K1),dpmdbet(nx-1,K1-1,IT,N), & dpjdbet(nx-1,K1,K1,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)),dht0(K1-1,K1,K1,K1), & temp0 do 2100 i=1,N if(nr(i).ge.2) then do 2102 t=2,nr(i) if(t.eq.2) then do 2104 jl=1,K1 temp0=0.d0 do 2106 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,1,il,t,i)*yv(jl)+ & gamz(2,1,il,t,i)*yv(jl)**2.d0) 2106 continue do 2108 ik=1,K1 if(ik.lt.K1) then h(1,jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,1,ik,t,i)*yv(jl)+ & gamz(2,1,ik,t,i)*yv(jl)**2)/(1.d0+temp0) else h(1,jl,K1,t,i)=1.d0/(1.d0+temp0) endif 2108 continue 2104 continue else do 2110 j1=1,K1 do 2110 j2=1,K1 temp0=0.d0 do 2112 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+ & gamz(1,1,il,t,i)*yv(j1)+gamz(2,1,il,t,i)*yv(j1)**2.d0+ & gamz(1,2,il,t,i)*yv(j2)+gamz(2,2,il,t,i)*yv(j2)**2.d0) 2112 continue do 2114 ik=1,K1 if(ik.lt.K1) then h(j2,j1,ik,t,i)=dexp(triangle(ik,t,i)+ & 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+temp0) else h(j2,j1,K1,t,i)=1.d0/(1.d0+temp0) endif 2114 continue 2110 continue endif 2102 continue endif 2100 continue do 2120 i=1,N if(nr(i).ge.2) then do 2122 t=2,nr(i) if(t.eq.2) then do 2124 k=1,K1-1 do 2124 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2125 j=1,K1 if(kg.eq.k) then dht0(kg,1,j,k)=h(1,j,k,t,i)*(1.d0-h(1,j,k,t,i)) else dht0(kg,1,j,k)=-h(1,j,k,t,i)*h(1,j,kg,t,i) endif dht0(kg,1,j,K1)=-h(1,j,K1,t,i)*h(1,j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,1,j,k)*pi(j,t-1,i) 2125 continue 2124 continue do 2126 k=1,K1-1 do 2126 ig=1,nx-1 psum(ig,k)=0.d0 2126 continue do 2128 k=1,K1-1 do 2130 j=1,K1 do 2130 ig=1,nx-1 if(j.eq.1) psum(ig,k)=psum(ig,k)+ & h(1,j,k,t,i)*dpmdbet(ig,j,t-1,i) if((j.gt.1).and.(j.lt.K1)) psum(ig,k)=psum(ig,k)+ & h(1,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)+ & h(1,j,k,t,i)*(-dpmdbet(ig,j-1,t-1,i)) 2130 continue do 2131 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) 2131 continue 2128 continue do 2132 ik=1,(nx-1)*(K1-1) do 2132 ig=1,(nx-1)*(K1-1) dhtsum1(ik,ig)=0.d0 2132 continue do 2133 k=1,K1-1 do 2133 kg=1,K1-1 do 2133 ig=1,nx-1 dhtsum1((k-1)*(nx-1)+ig,(kg-1)*(nx-1)+ig)=dhtsum(k,kg) 2133 continue do 2134 ik=1,K1-1 do 2134 l=1,nx-1 psum1((ik-1)*(nx-1)+l)=psum(l,ik) 2134 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 2135 k=1,K1-1 do 2135 ig=1,nx-1 dtrdbet(ig,k,t,i)=res((k-1)*(nx-1)+ig) 2135 continue do 2137 j1=1,K1 do 2137 j2=1,K1 do 2137 ig=1,nx-1 dpjdbet(ig,j2,j1,t,i)=0.d0 do 2138 l=1,K1-1 dpjdbet(ig,j2,j1,t,i)=dpjdbet(ig,j2,j1,t,i)+ & dht0(l,1,j2,j1)*dtrdbet(ig,l,t,i)*pi(j2,t-1,i) 2138 continue if(j2.eq.1) then dpjdbet(ig,j2,j1,t,i)=dpjdbet(ig,j2,j1,t,i)+ & h(1,j2,j1,t,i)*dpmdbet(ig,j2,t-1,i) else if(j2.eq.K1) then dpjdbet(ig,j2,j1,t,i)=dpjdbet(ig,j2,j1,t,i)- & h(1,j2,j1,t,i)*dpmdbet(ig,j2-1,t-1,i) else dpjdbet(ig,j2,j1,t,i)=dpjdbet(ig,j2,j1,t,i)+ & h(1,j2,j1,t,i)* & (dpmdbet(ig,j2,t-1,i)-dpmdbet(ig,j2-1,t-1,i)) endif 2137 continue else do 2140 k=1,K1-1 do 2140 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2142 j1=1,K1 do 2142 j2=1,K1 if(kg.eq.k) then dht0(kg,j2,j1,k)=h(j2,j1,k,t,i)*(1.d0-h(j2,j1,k,t,i)) else dht0(kg,j2,j1,k)=-h(j2,j1,k,t,i)*h(j2,j1,kg,t,i) endif dht0(kg,j2,j1,K1)=-h(j2,j1,K1,t,i)*h(j2,j1,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,j2,j1,k)*pj(j2,j1,t-1,i) 2142 continue 2140 continue do 2144 k=1,K1-1 do 2144 ig=1,nx-1 psum(ig,k)=0.d0 2144 continue do 2146 k=1,K1-1 do 2148 ig=1,nx-1 do 2148 j1=1,K1 do 2148 j2=1,K1 psum(ig,k)=psum(ig,k)+ & h(j2,j1,k,t,i)*dpjdbet(ig,j2,j1,t-1,i) 2148 continue do 2149 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) 2149 continue 2146 continue do 2150 ik=1,(nx-1)*(K1-1) do 2150 ig=1,(nx-1)*(K1-1) dhtsum1(ik,ig)=0.d0 2150 continue do 2152 k=1,K1-1 do 2152 kg=1,K1-1 do 2152 ig=1,nx-1 dhtsum1((k-1)*(nx-1)+ig,(kg-1)*(nx-1)+ig)=dhtsum(k,kg) 2152 continue do 2154 ik=1,K1-1 do 2154 l=1,nx-1 psum1((ik-1)*(nx-1)+l)=psum(l,ik) 2154 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 2156 k=1,K1-1 do 2156 ig=1,nx-1 dtrdbet(ig,k,t,i)=res((k-1)*(nx-1)+ig) 2156 continue do 2157 j1=1,K1 do 2157 j2=1,K1 do 2157 ig=1,nx-1 dpjdbet(ig,j2,j1,t,i)=0.d0 do 2158 j3=1,K1 do 2159 l=1,K1-1 dpjdbet(ig,j2,j1,t,i)=dpjdbet(ig,j2,j1,t,i)+ & dht0(l,j3,j2,j1)*dtrdbet(ig,l,t,i)*pj(j3,j2,t-1,i) 2159 continue dpjdbet(ig,j2,j1,t,i)=dpjdbet(ig,j2,j1,t,i)+ & h(j3,j2,j1,t,i)*dpjdbet(ig,j3,j2,t-1,i) 2158 continue 2157 continue endif 2122 continue endif 2120 continue return end c Calculate d tri / d beta_0 subroutine dtrbet0(gamz,triangle,yv,pi,pj,dpmdbet0, & N,IT,nr,K1,norder,ialph,dtrdbet0) integer t,nr(N) double precision pj(K1,K1,IT,N), & h(K1,K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,norder,K1-1,IT,N), & triangle(K1-1,IT,N),yv(K1),dpmdbet0(K1-1,K1-1,IT,N), & dpjdbet0(K1-1,K1,K1,IT,N), & dhtsum(K1-1,K1-1),dtrdbet0(K1-1,K1-1,IT,N), & psum(K1-1),Ainv(K1-1,K1-1),resbet0(K1-1), & dht0(K1-1,K1,K1,K1), & temp0 do 2200 i=1,N if(nr(i).ge.2) then do 2202 t=2,nr(i) if(t.eq.2) then do 2204 jl=1,K1 temp0=0.d0 do 2206 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,1,il,t,i)*yv(jl)+ & gamz(2,1,il,t,i)*yv(jl)**2.d0) 2206 continue do 2208 ik=1,K1 if(ik.lt.K1) then h(1,jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,1,ik,t,i)*yv(jl)+ & gamz(2,1,ik,t,i)*yv(jl)**2)/(1.d0+temp0) else h(1,jl,K1,t,i)=1.d0/(1.d0+temp0) endif 2208 continue 2204 continue else do 2210 j1=1,K1 do 2210 j2=1,K1 temp0=0.d0 do 2212 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+ & gamz(1,1,il,t,i)*yv(j1)+gamz(2,1,il,t,i)*yv(j1)**2.d0+ & gamz(1,2,il,t,i)*yv(j2)+gamz(2,2,il,t,i)*yv(j2)**2.d0) 2212 continue do 2214 ik=1,K1 if(ik.lt.K1) then h(j2,j1,ik,t,i)=dexp(triangle(ik,t,i)+ & 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+temp0) else h(j2,j1,K1,t,i)=1.d0/(1.d0+temp0) endif 2214 continue 2210 continue endif 2202 continue endif 2200 continue do 2220 i=1,N if(nr(i).ge.2) then do 2222 t=2,nr(i) if(t.eq.2) then do 2224 k=1,K1-1 do 2224 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2225 j=1,K1 if(kg.eq.k) then dht0(kg,1,j,k)=h(1,j,k,t,i)*(1.d0-h(1,j,k,t,i)) else dht0(kg,1,j,k)=-h(1,j,k,t,i)*h(1,j,kg,t,i) endif dht0(kg,1,j,K1)=-h(1,j,K1,t,i)*h(1,j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,1,j,k)*pi(j,t-1,i) 2225 continue 2224 continue do 2226 ia=1,K1-1 do 2228 k=1,K1-1 psum(k)=0.d0 do 2230 j=1,K1 if(j.eq.1) psum(k)=psum(k)+ & h(1,j,k,t,i)*dpmdbet0(ia,j,t-1,i) if((j.gt.1).and.(j.lt.K1)) psum(k)=psum(k)+ & h(1,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)+ & h(1,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) c if(k.eq.K1) psum(k)=-dpmdbet0(ia,k-1,t,i)-psum(k) 2228 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 2232 ik=1,K1-1 dtrdbet0(ia,ik,t,i)=resbet0(ik) 2232 continue 2226 continue do 2233 j1=1,K1 do 2233 j2=1,K1 do 2233 ig=1,K1-1 dpjdbet0(ig,j2,j1,t,i)=0.d0 do 2234 l=1,K1-1 dpjdbet0(ig,j2,j1,t,i)=dpjdbet0(ig,j2,j1,t,i)+ & dht0(l,1,j2,j1)*dtrdbet0(ig,l,t,i)*pi(j2,t-1,i) 2234 continue if(j2.eq.1) then dpjdbet0(ig,j2,j1,t,i)=dpjdbet0(ig,j2,j1,t,i)+ & h(1,j2,j1,t,i)*dpmdbet0(ig,j2,t-1,i) else if(j2.eq.K1) then dpjdbet0(ig,j2,j1,t,i)=dpjdbet0(ig,j2,j1,t,i)- & h(1,j2,j1,t,i)*dpmdbet0(ig,j2-1,t-1,i) else dpjdbet0(ig,j2,j1,t,i)=dpjdbet0(ig,j2,j1,t,i)+ & h(1,j2,j1,t,i)* & (dpmdbet0(ig,j2,t-1,i)-dpmdbet0(ig,j2-1,t-1,i)) endif 2233 continue else do 2240 k=1,K1-1 do 2240 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2242 j1=1,K1 do 2242 j2=1,K1 if(kg.eq.k) then dht0(kg,j2,j1,k)=h(j2,j1,k,t,i)*(1.d0-h(j2,j1,k,t,i)) else dht0(kg,j2,j1,k)=-h(j2,j1,k,t,i)*h(j2,j1,kg,t,i) endif dht0(kg,j2,j1,K1)=-h(j2,j1,K1,t,i)*h(j2,j1,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,j2,j1,k)*pj(j2,j1,t-1,i) 2242 continue 2240 continue do 2246 ia=1,K1-1 do 2248 k=1,K1-1 psum(k)=0.d0 do 2250 j1=1,K1 do 2250 j2=1,K1 psum(k)=psum(k)+ & h(j2,j1,k,t,i)*dpjdbet0(ia,j2,j1,t-1,i) 2250 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) c if(k.eq.K1) psum(k)=-dpmdbet0(ia,k-1,t,i)-psum(k) 2248 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 2252 ik=1,K1-1 dtrdbet0(ia,ik,t,i)=resbet0(ik) 2252 continue 2246 continue do 2257 j1=1,K1 do 2257 j2=1,K1 do 2257 ig=1,K1-1 dpjdbet0(ig,j2,j1,t,i)=0.d0 do 2258 j3=1,K1 do 2259 l=1,K1-1 dpjdbet0(ig,j2,j1,t,i)=dpjdbet0(ig,j2,j1,t,i)+ & dht0(l,j3,j2,j1)*dtrdbet0(ig,l,t,i)*pj(j3,j2,t-1,i) 2259 continue dpjdbet0(ig,j2,j1,t,i)=dpjdbet0(ig,j2,j1,t,i)+ & h(j3,j2,j1,t,i)*dpjdbet0(ig,j3,j2,t-1,i) 2258 continue 2257 continue endif 2222 continue endif 2220 continue return end c Calculate d tri / d alpha_1b^(a) subroutine dtralp1(z,triangle,yv,pi,pj,dhdgam1,gamz, & N,IT,nr,K1,nalpha,ialph,norder,dtrdalp1) integer t,nr(N) double precision z(nalpha,IT,N),yv(K1), & pj(K1,K1,IT,N), & h(K1,K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,norder,K1-1,IT,N), & triangle(K1-1,IT,N), & dhdgam1(ialph,K1-1,K1,K1,K1,IT,N), & dhtsum(K1-1,K1-1),dtrdalp1(nalpha,ialph,K1-1,K1-1,IT,N), & psum(nalpha,K1-1),dht0(K1-1,K1,K1,K1), & Ainv((K1-1)*nalpha,(K1-1)*nalpha),res((K1-1)*nalpha), & dhtsum1(nalpha*(K1-1),nalpha*(K1-1)), & dpjdalp1(nalpha,ialph,K1-1,K1,K1,IT,N), & psum1((K1-1)*nalpha),temp0 do 2299 i=1,N do 2299 t=1,IT do 2299 k=1,K1-1 do 2299 ia=1,K1-1 do 2299 ib=1,ialph do 2299 ig=1,nalpha dtrdalp1(ig,ib,ia,k,t,i)=0.d0 2299 continue do 2300 i=1,N if(nr(i).ge.2) then do 2302 t=2,nr(i) if(t.eq.2) then do 2304 jl=1,K1 temp0=0.d0 do 2306 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+gamz(1,1,il,t,i)*yv(jl)+ & gamz(2,1,il,t,i)*yv(jl)**2.d0) 2306 continue do 2308 ik=1,K1 if(ik.lt.K1) then h(1,jl,ik,t,i)=dexp(triangle(ik,t,i)+ & gamz(1,1,ik,t,i)*yv(jl)+ & gamz(2,1,ik,t,i)*yv(jl)**2.d0)/(1.d0+temp0) else h(1,jl,K1,t,i)=1.d0/(1.d0+temp0) endif 2308 continue 2304 continue else do 2310 j1=1,K1 do 2310 j2=1,K1 temp0=0.d0 do 2312 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+ & gamz(1,1,il,t,i)*yv(j1)+gamz(2,1,il,t,i)*yv(j1)**2.d0+ & gamz(1,2,il,t,i)*yv(j2)+gamz(2,2,il,t,i)*yv(j2)**2.d0) 2312 continue do 2314 ik=1,K1 if(ik.lt.K1) then h(j2,j1,ik,t,i)=dexp(triangle(ik,t,i)+ & 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+temp0) else h(j2,j1,K1,t,i)=1.d0/(1.d0+temp0) endif 2314 continue 2310 continue endif 2302 continue endif 2300 continue do 2320 i=1,N if(nr(i).ge.2) then do 2322 t=2,nr(i) if(t.eq.2) then do 2324 k=1,K1-1 do 2324 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2326 j=1,K1 dht0(kg,2,j,k)=0.d0 dht0(kg,2,j,K1)=0.d0 if(kg.eq.k) dht0(kg,1,j,k)= & h(1,j,k,t,i)*(1.d0-h(1,j,k,t,i)) if(kg.ne.k) dht0(kg,1,j,k)=-h(1,j,k,t,i)*h(1,j,kg,t,i) dht0(kg,1,j,K1)=-h(1,j,K1,t,i)*h(1,j,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,1,j,k)*pi(j,t-1,i) 2326 continue 2324 continue do 2328 ia=1,K1-1 do 2328 ib=1,ialph do 2332 k=1,K1-1 do 2332 in=1,nalpha psum(in,k)=0.d0 do 2334 ig=1,K1 psum(in,k)=psum(in,k)- & dhdgam1(ib,ia,1,ig,k,t,i)*z(in,t,i)*pi(ig,t-1,i) 2334 continue 2332 continue do 2336 ik=1,nalpha*(K1-1) do 2336 ig=1,nalpha*(K1-1) dhtsum1(ik,ig)=0.d0 2336 continue do 2338 k=1,K1-1 do 2338 kg=1,K1-1 do 2338 ig=1,nalpha dhtsum1((k-1)*nalpha+ig,(kg-1)*nalpha+ig)=dhtsum(k,kg) 2338 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 2342 k=1,K1-1 do 2344 ig=1,nalpha dtrdalp1(ig,ib,ia,k,t,i)=res((k-1)*nalpha+ig) 2344 continue 2342 continue 2328 continue do 2345 j1=1,K1 do 2345 j2=1,K1 do 2345 ia=1,K1-1 do 2345 ib=1,ialph do 2345 ig=1,nalpha dpjdalp1(ig,ib,ia,j2,j1,t,i)=0.d0 do 2346 l=1,K1-1 dpjdalp1(ig,ib,ia,j2,j1,t,i)=dpjdalp1(ig,ib,ia,j2,j1,t,i)+ & dht0(l,1,j2,j1)*dtrdalp1(ig,ib,ia,l,t,i)*pi(j2,t-1,i) 2346 continue dpjdalp1(ig,ib,ia,j2,j1,t,i)=dpjdalp1(ig,ib,ia,j2,j1,t,i)+ & dhdgam1(ib,ia,1,j2,j1,t,i)*z(ig,t,i)*pi(j2,t-1,i) 2345 continue else do 2350 k=1,K1-1 do 2350 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2352 j1=1,K1 do 2352 j2=1,K1 if(kg.eq.k) dht0(kg,j2,j1,k)= & h(j2,j1,k,t,i)*(1.d0-h(j2,j1,k,t,i)) if(kg.ne.k) dht0(kg,j1,j2,k)= & -h(j2,j1,k,t,i)*h(j2,j1,kg,t,i) dht0(kg,j2,j1,K1)=-h(j2,j1,K1,t,i)*h(j2,j1,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,j2,j1,k)*pj(j2,j1,t-1,i) 2352 continue 2350 continue do 2358 ia=1,K1-1 do 2358 ib=1,ialph do 2362 k=1,K1-1 do 2362 in=1,nalpha psum(in,k)=0.d0 do 2364 j1=1,K1 do 2364 j2=1,K1 psum(in,k)=psum(in,k)- & dhdgam1(ib,ia,j2,j1,k,t,i)*z(in,t,i)*pj(j2,j1,t-1,i)- & h(j2,j1,k,t,i)*dpjdalp1(in,ib,ia,j2,j1,t-1,i) 2364 continue 2362 continue do 2366 ik=1,nalpha*(K1-1) do 2366 ig=1,nalpha*(K1-1) dhtsum1(ik,ig)=0.d0 2366 continue do 2368 k=1,K1-1 do 2368 kg=1,K1-1 do 2368 ig=1,nalpha dhtsum1((k-1)*nalpha+ig,(kg-1)*nalpha+ig)=dhtsum(k,kg) 2368 continue do 2370 ik=1,K1-1 do 2370 l=1,nalpha psum1((ik-1)*nalpha+l)=psum(l,ik) 2370 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 2372 k=1,K1-1 do 2374 ig=1,nalpha dtrdalp1(ig,ib,ia,k,t,i)=res((k-1)*nalpha+ig) 2374 continue 2372 continue 2358 continue do 2380 j1=1,K1 do 2380 j2=1,K1 do 2380 ia=1,K1-1 do 2380 ib=1,ialph do 2380 ig=1,nalpha dpjdalp1(ig,ib,ia,j2,j1,t,i)=0.d0 do 2381 j3=1,K1 do 2382 l=1,K1-1 dpjdalp1(ig,ib,ia,j2,j1,t,i)=dpjdalp1(ig,ib,ia,j2,j1,t,i)+ & dht0(l,j3,j2,j1)*dtrdalp1(ig,ib,ia,l,t,i)*pj(j3,j2,t-1,i) 2382 continue dpjdalp1(ig,ib,ia,j2,j1,t,i)=dpjdalp1(ig,ib,ia,j2,j1,t,i)+ & dhdgam1(ib,ia,j3,j2,j1,t,i)*z(ig,t,i)*pj(j3,j2,t-1,i)+ & h(j3,j2,j1,t,i)*dpjdalp1(ig,ib,ia,j3,j2,t-1,i) 2381 continue 2380 continue endif 2322 continue endif 2320 continue return end c Calculate d tri / d alpha_2b^(a) subroutine dtralp2(z,triangle,yv,pi,pj,dhdgam2,gamz, & N,IT,nr,K1,nalpha,ialph,norder,dtrdalp2) integer t,nr(N) double precision z(nalpha,IT,N),yv(K1), & pj(K1,K1,IT,N), & h(K1,K1,K1,IT,N),pi(K1,IT,N),gamz(ialph,norder,K1-1,IT,N), & triangle(K1-1,IT,N), & dhdgam2(ialph,K1-1,K1,K1,K1,IT,N), & dhtsum(K1-1,K1-1),dtrdalp2(nalpha,ialph,K1-1,K1-1,IT,N), & psum(nalpha,K1-1),dht0(K1-1,K1,K1,K1), & Ainv((K1-1)*nalpha,(K1-1)*nalpha),res((K1-1)*nalpha), & dhtsum1(nalpha*(K1-1),nalpha*(K1-1)), & psum1((K1-1)*nalpha), & dpjdalp2(nalpha,ialph,K1-1,K1,K1,IT,N),temp0 do 2399 i=1,N do 2399 t=1,IT do 2399 k=1,K1-1 do 2399 ia=1,K1-1 do 2399 ib=1,ialph do 2399 ig=1,nalpha dtrdalp2(ig,ib,ia,k,t,i)=0.d0 2399 continue do 2400 i=1,N if(nr(i).ge.3) then do 2402 t=3,nr(i) do 2410 j1=1,K1 do 2410 j2=1,K1 temp0=0.d0 do 2412 il=1,K1-1 temp0=temp0+dexp(triangle(il,t,i)+ & gamz(1,1,il,t,i)*yv(j1)+gamz(2,1,il,t,i)*yv(j1)**2.d0+ & gamz(1,2,il,t,i)*yv(j2)+gamz(2,2,il,t,i)*yv(j2)**2.d0) 2412 continue do 2414 ik=1,K1 if(ik.lt.K1) then h(j2,j1,ik,t,i)=dexp(triangle(ik,t,i)+ & 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+temp0) else h(j2,j1,K1,t,i)=1.d0/(1.d0+temp0) endif 2414 continue 2410 continue 2402 continue endif 2400 continue do 2420 i=1,N if(nr(i).ge.2) then do 2422 t=3,nr(i) do 2442 k=1,K1-1 do 2442 ia=1,K1-1 do 2442 ib=1,ialph do 2442 ig=1,nalpha dtrdalp2(ig,ib,ia,k,2,i)=0.d0 2442 continue do 2445 j1=1,K1 do 2445 j2=1,K1 do 2445 ia=1,K1-1 do 2445 ib=1,norder do 2445 ig=1,nalpha dpjdalp2(ig,ib,ia,j2,j1,2,i)=0.d0 2445 continue do 2450 k=1,K1-1 do 2450 kg=1,K1-1 dhtsum(k,kg)=0.d0 do 2452 j1=1,K1 do 2452 j2=1,K1 if(kg.eq.k) dht0(kg,j2,j1,k)= & h(j2,j1,k,t,i)*(1.d0-h(j2,j1,k,t,i)) if(kg.ne.k) dht0(kg,j1,j2,k)= & -h(j2,j1,k,t,i)*h(j2,j1,kg,t,i) dht0(kg,j2,j1,K1)=-h(j2,j1,K1,t,i)*h(j2,j1,kg,t,i) dhtsum(k,kg)=dhtsum(k,kg)+dht0(kg,j2,j1,k)*pj(j2,j1,t-1,i) 2452 continue 2450 continue do 2458 ia=1,K1-1 do 2458 ib=1,ialph do 2460 k=1,K1-1 do 2460 in=1,nalpha psum(in,k)=0.d0 do 2461 j1=1,K1 do 2461 j2=1,K1 psum(in,k)=psum(in,k)- & dhdgam2(ib,ia,j2,j1,k,t,i)*z(in,t,i)*pj(j2,j1,t-1,i)- & h(j2,j1,k,t,i)*dpjdalp2(in,ib,ia,j2,j1,t-1,i) 2461 continue 2460 continue do 2466 ik=1,nalpha*(K1-1) do 2466 ig=1,nalpha*(K1-1) dhtsum1(ik,ig)=0.d0 2466 continue do 2468 k=1,K1-1 do 2468 kg=1,K1-1 do 2468 ig=1,nalpha dhtsum1((k-1)*nalpha+ig,(kg-1)*nalpha+ig)=dhtsum(k,kg) 2468 continue do 2470 ik=1,K1-1 do 2470 l=1,nalpha psum1((ik-1)*nalpha+l)=psum(l,ik) 2470 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 2472 k=1,K1-1 do 2474 ig=1,nalpha dtrdalp2(ig,ib,ia,k,t,i)=res((k-1)*nalpha+ig) 2474 continue 2472 continue 2458 continue do 2480 j1=1,K1 do 2480 j2=1,K1 do 2480 ia=1,K1-1 do 2480 ib=1,ialph do 2480 ig=1,nalpha dpjdalp2(ig,ib,ia,j2,j1,t,i)=0.d0 do 2481 j3=1,K1 do 2482 l=1,K1-1 dpjdalp2(ig,ib,ia,j2,j1,t,i)=dpjdalp2(ig,ib,ia,j2,j1,t,i)+ & dht0(l,j3,j2,j1)*dtrdalp2(ig,ib,ia,l,t,i)*pj(j3,j2,t-1,i) 2482 continue dpjdalp2(ig,ib,ia,j2,j1,t,i)=dpjdalp2(ig,ib,ia,j2,j1,t,i)+ & dhdgam2(ib,ia,j3,j2,j1,t,i)*z(ig,t,i)*pj(j3,j2,t-1,i)+ & h(j3,j2,j1,t,i)*dpjdalp2(ig,ib,ia,j3,j2,t-1,i) 2481 continue 2480 continue 2422 continue endif 2420 continue return end