Biomedical Engineering Reference
In-Depth Information
p[i,]=y$p;Tnt[i,]=y$Tnt;Tt[i,]=y$Tt
}
#####Step:E-step2:E[R]E[Rsq]E[expR]##################
for(iin1:N)
{
for(jin1:GH_n)
{
temp=sqrt(2*sigma2)*Xj[j]
L[j]=exp(sum(p[i,]*log(al01)+(dij[i,]*d[i])*log(al02)-
al01*Tnt[i,]*exp(beta1*z[i]+temp)-al02*Tnt[i,]*
exp(beta2*z[i]+temp)-al02*Tt[i,]*exp(beta3*z[i]+temp))+
(de[i]*beta1+(1-de[i])*d[i]*beta2+
de[i]*d[i]*beta3)*z[i]+(de[i]+d[i])*temp)
f0[j]=Wj[j]*(1/sqrt(pi))*L[j]#
f1[j]=Wj[j]*(1/sqrt(pi))*L[j]*temp#:E(R)
f2[j]=Wj[j]*(1/sqrt(pi))*L[j]*(exp(temp))#:E(exp(R))
f3[j]=Wj[j]*(1/sqrt(pi))*L[j]*(temp)^2#:E(R^2)
}
if(sum(which(f1== ' NaN ' ))==0&sum(which(f1== ' Inf ' ))==0)
{
R[i]=sum(f1)/sum(f0)
expR[i]=sum(f2)/sum(f0)
Rsq[i]=sum(f3)/sum(f0)
}else
{flag=0;sim_flag=1;break}
}
#####Step:M-step:NewtonRaphson######################
res=c(beta1,beta2,beta3)#notusingmaxNRfunction
res=res-gr_beta(res)%*%(qr.solve(he_beta(res),tol=1e-17))#
beta1=res[1];beta2=res[2];beta3=res[3]#
al01=colSums(p)/(colSums(Tnt*exp(beta1*z)*expR))
al02=colSums(dij*d)/(colSums(Tnt*exp(beta2*z)*expR)+
colSums(Tt*exp(beta3*z)*expR))
 
Search WWH ::




Custom Search