先描述一下具体问题,首先是目标函数

(代码里的目标函数没有参数贝塔)其中,对r[t-x]的参数估计为

,a、k的参数估计公式类似
下面是我的代码:
#likelihood法参数估计
###(1)set the start value
a<-c(rep(0,90));
k<-c(rep(1,T));
r<-c(rep(1,89+T));
###(2)construct the functions
#####(2.1)construct the log-likelihood function
logL<-function(a,k,r){
L<-matrix(rep(NA,90*T),90,T);
for (x in 1:90){
for (t in 1:T){
L[x,t]<-D[x,t]*(a[x]+k[t]+r[t-x])-E[x,t]*exp(a[x]+k[t]+r[t-x]);
}
}
logL<-sum(L);
return(logL);
}
#####(2.2)the estimated number od deaths
dd<-function(a,k,r){
d<-matrix(rep(0,90*T),90,T);#the estimated number of deaths in each age group and each year
for(t in 1:T){
for(x in 1:90){
d[x,t]<-E[x,t]*exp(a[x]+k[t]+r[t-x]);
}
}
return(d);
}
#####(2.3)set the iteration rule of a(x)
aa<-function(a,k,r){
d<-dd(a,k,r);
for(x in 1:90){
a[x]<-a[x]-sum(D[x,]-d[x,])/(-sum(d[x,]));
}
return(a);
}
#####(2.4)set the iteration rule of k(t)
kk<-function(a,k,r){
d<-dd(a,k,r);
for(t in 1:T){
k[t]<-k[t]+sum(D[,t]-d[,t])/(sum(d[,t]));
}
return(k);
}
#####(2.5)set the iteration rule of r(z)
rr<-function(a,k,r){
d<-dd(a,k,r);
for (t in 1:T)
for (x in 1:90){
r[t-x]<-r[t-x]+sum((D[x,t]-d[x,t]))/(sum(d[x,t]));
}
return(r);
}
###(3)iteration process
j=1;
diff_L<-logL(a,k,r);
while(abs(diff_L)>10^(-10)){
logL0<-logL(a,b,k);
if(j%%3==1){a<-aa(a,k,r);}
if(j%%3==2){k<-kk(a,k,r);}
if(j%%3==0){r<-rr(a,k,r);}
logL1<-logL(a,k,r);
diff_L<-logL1-logL0;
j=j+1;
}
想问一下各位大佬,r[t-x]的参数估计应该用R语言怎么写(t、x都是变量,但是只要t-x相等,r[t-x]就相等)
还有就是我在循环里的diff_L<-logL(a,k,r); 这一句它总报错,说是replacement has length zero