zonghedata1<-zonghedata
break_year<-2010
T<-break_year-1970+1
n<-2014-break_year
D_all<-matrix(rep(0,17*41), 17,41)
E_all<-matrix(rep(0,17*41),17,41)
m_all<-matrix(rep(0,17*41),17,41)
lnm_all<-matrix(rep(0,17*41),17,41)
for(t in 1:T){
  for (x in 1:17)
  {D_all[x,t]<-zonghedata1[x+17*(t-1),3]
  E_all[x,t]<-zonghedata1[x+17*(t-1),4]
  m_all[x,t]<-D_all[x,t]/E_all[x,t]
  lnm_all[x,t]<-log(m_all[x,t])
  }
}
D<-D_all[,1:t]
E<-E_all[,1:t]
m<-m_all[,1:t]
lnm<-lnm_all[,1:t]
age<-zonghedata1[1:17,2]
year<-as.numeric(levels(factor(zonghedata1[,1])))
#step1:set the start value for a0,b0,k0
a<-c(rep(0,17))
b<-c(rep(1,17))
k<-c(rep(0,41))
#step2: construct the log_likelihood function
logL<-function(a,b,k){
  L<-matrix(rep(NA,17*t),17,t)
  for (x in 1:17) {for(t in 1:T){
    L[x,t]<-D[x,t]*log(1-exp(-exp(a[x]+b[x]*k[t])))-D[x,t]*exp(a[x]+b[x]*k[t])####对数似然函数
  }}
  logL<-sum(L)
  return(logL)
}
#step3:the estimated number of deaths
dd<-function(a,b,k){
  q<-matrix(rep(0,17*t),17,t)
  for(t in 1:T){
    for (x in 1:17) {
      q[x,t]<-1-exp(-exp(a[x]+b[x]*k[t]))
    }
  }
  return(q)
}
#step4:set the iteration rule of a(x)
aa<-function(a,b,k){
  q<-dd(a,b,k)
  for (x in 1:17){
  a[x]<-a[x]-sum((exp(a[x]-b[x]*k[t]))*(D[x,]/q[x,]-E[x,]))/sum((exp(a[x]-b[x]*k[t])*((D[x,]*q[x,]+D[x,]*exp(a[x]+b[x]*k[t])*exp(-exp(a[x]+b[x]*k[t]))))/(q[x,])^2-E[x,]))
  }
  return(a)
}
# step5:set the iteration rule of k(t)
kk<-function(a,b,k){
  q<-dd(a,b,k)
  for (t in 1:T) {k[t]<-k[t]-sum((b*exp(a[x]+b[x]*k[t]))*(D[,t]/q[,t]-E[,t]))/sum((b^2)*exp(a[x]+b[x]*k[t])*(((D[,t]*q[,t]+D[,t]*exp(a[x]+b[x]*k[t])*exp(-exp(a[x]+b[x]*k[t]))))/(q[,t])^2-E[,t]))
  }
  return(k)
}
#step6: set the iteration of b(x)
bb<- function(a,b,k){
  q<-dd(a,b,k)
  for (x in 1:17) {b[x]<-b[x]-sum((k*exp(a[x]+b[x]*k[t]))*(D[x,]/q[x,]-E[x,]))/sum((k^2)*exp(a[x]+b[x]*k[t])*(((D[x,]*q[x,]+D[x,]*exp(a[x]+b[x]*k[t])*exp(-exp(a[x]+b[x]*k[t]))))/(q[x,])^2-E[x,]))
  }
  return(b)
}
#step7:iteration process
j=1
diff_L<-logL(a,b,k)
while(abs(diff_L)>10^(-10))
{logL0<-logL(a,b,k)
if(j%%3==1){a<-aa(a,b,k)}
if(j%%3==2){k<-kk(a,b,k)}
if(j%%3==0){b<-bb(a,b,k)}
logL1<-logL(a,b,k)
diff_L<-logL1-logL0
j=j+1
}
Error in while (abs(diff_L) > 10^(-10)) { :
  missing value where TRUE/FALSE needed
#step8:parameter adjustment
a1<-a
b1<-b
k1<-k
c1<-sum(b1)
c2<-c1*sum(k1)/T
A_MLE<-a1+c2*b/c1
B_MLE<-b1/c1
K_MLE<-c1*k1-c2


提示我出现缺失值，找了好久不都不知道为什么