diff --git a/R/kfino.R b/R/kfino.R index 87678b30fefe2524c71e01392d2655b3d1209965..2291d883ae2197065f4a705d3a88f47957911dfa 100644 --- a/R/kfino.R +++ b/R/kfino.R @@ -287,32 +287,70 @@ kfino_fit<-function(datain,Tvar,Yvar, bornem0=quantile(Y[1:N/4], probs = c(.2, .8)) m0opt=quantile(Y[1:N/4], probs = c(.5)) mmopt=quantile(Y[(3*N/4):N], probs = c(.5)) - cat("range m0: ",bornem0,"\n") cat("initial m0opt: ",m0opt,"\n") cat("initial mmopt: ",mmopt,"\n") - popt=0.5 - EMopt=KBO_EM(param=list(m0=m0opt, - mm=mmopt, - pp=popt, - aa=aa, - expertMin=expertMin, - expertMax=expertMax, - sigma2_mm=sigma2_mm, - sigma2_m0=sigma2_m0, - sigma2_pp=sigma2_pp, - K=K), - Y=Y,Tps=Tps,N=N,dix=dix,kappaOpt=kappaOpt) + # 1er tour + Vopt=KBO_EM(param=list(m0=m0opt, + mm=mmopt, + pp=popt, + aa=aa, + expertMin=expertMin, + expertMax=expertMax, + sigma2_mm=sigma2_mm, + sigma2_m0=sigma2_m0, + sigma2_pp=sigma2_pp, + K=K), + Y=Y,Tps=Tps,N=N,dix=dix,kappaOpt=kappaOpt)$likelihood + + # boucle + N_etape_EM=10 + t=1 + for (m0_tmp in seq(bornem0[1],bornem0[2],2) ){ + print(t) + t=t+1 + m_tmp=m0_tmp + p_tmp=0.5 + for (k in 1:N_etape_EM){ + Res_EM=KBO_EM(param=list(mm=m_tmp, + pp=p_tmp, + m0=m0_tmp, + aa=aa, + expertMin=expertMin, + expertMax=expertMax, + sigma2_m0=sigma2_m0, + sigma2_mm=sigma2_mm, + sigma2_pp=sigma2_pp, + K=K), + kappaOpt=kappaOpt, Y=Y,Tps=Tps,N=N,dix=dix) + m0_tmp=Res_EM$m0[[1]] + m_tmp=Res_EM$mm[[1]] + p_tmp=Res_EM$pp + } + V=Res_EM$likelihood + + #print(V) + print(m0_tmp) + print(m_tmp) + if (V > Vopt){ + print("EM converged.") + Vopt=V + m0opt=m0_tmp + mmopt=m_tmp + popt=p_tmp + } + } + print("Optimised parameters with EM method: ") - cat("Optimized m0: ",EMopt$m0[[1]],"\n") - cat("Optimized mm: ",EMopt$mm[[1]],"\n") - cat("Optimized pp: ",EMopt$pp[[1]],"\n") + cat("Optimized m0: ",m0opt,"\n") + cat("Optimized mm: ",mmopt,"\n") + cat("Optimized pp: ",popt,"\n") print("-------:") - resultat=KBO_known(param=list(mm=EMopt$mm[[1]], - pp=EMopt$pp[[1]], - m0=EMopt$m0[[1]], + resultat=KBO_known(param=list(mm=mmopt, + pp=popt, + m0=m0opt, aa=aa, expertMin=expertMin, expertMax=expertMax,