Biomedical Engineering Reference
In-Depth Information
convLik=matrix(NA,nrow=nMax,ncol=1)
#E-M################################################################
w=0;flag=0;#w:countofEMiteration,flag:stopindicator
error_flag=0;#error_flag:errorindicator
#cat("LLbeta1beta2beta3theta","\n")
while(w<nMax&flag==0)
{
#####previousdatasaving#############################################
pal01=al01;pal02=al02;pbeta1=beta1;
pbeta2=beta2;pbeta3=beta3;ptheta=theta
#####Step:E-step1:E[Nt]E[Nnt]process#############################
for(iin1:length(t))
{
y=q_s(t[i],dij[i,],de[i],d[i],interval,al01*
exp(beta1*z[i])*R[i],al02*exp(beta2*z[i])*R[i],al02*exp(beta3*z[i])*R[i])
p[i,]=y$p;Tnt[i,]=y$Tnt;Tt[i,]=y$Tt
}
#####Step:E-step2:E[R]E[logR]#####################################
A=1/theta+de+d;
C=1/theta+rowSums((Tnt*exp(beta1*z)*R)%*%al01)+
rowSums((Tnt*exp(beta2*z)*R)%*%al02)+rowSums((Tt*exp(beta3*z)*R)%*%al02)
R=A/C
logR=digamma(A)-log(C)
#####Step:M-step:NewtonRaphson####################################
res=c(beta1,beta2,beta3)#notusingmaxNRfunction
res=res-gr_beta(res)%*%(qr.solve(he_beta(res),tol=1e-30));#
beta1=res[1];beta2=res[2];beta3=res[3]#
al01=colSums(p)/(colSums(Tnt*exp(beta1*z)*R))
al02=colSums(dij*d)/(colSums(Tnt*exp(beta2*z)*R)
+colSums(Tt*exp(beta3*z)*R))
theta=theta-gr_theta(theta)/he_theta(theta)
if(w==nMax){flag=1;error_flag=1}else#error
{
Search WWH ::
Custom Search