MaJ execution de KBO_EM() dans kfino_fit() - par dichotomie

R/kfino.R

@@ -279,42 +279,25 @@ kfino_fit<-function(datain,Tvar,Yvar,
                          threshold=threshold,Y=Y,Tps=Tps,N=N,kappa=kappa)
 
     } else if (N > 50){
-      # optimization without sub-sampling
+      # optimization without sub-sampling, 2 methods, EM or ML
       if (method == "EM"){
-        print("-------:")
-        print("Optimisation of initial parameters with EM method - result:")
-        print("no sub-sampling performed:")
-        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")
-        popt=0.5
-
-        # 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) ){
-          t=t+1
+          print("-------:")
+          print("Optimisation of initial parameters with EM method - result:")
+          print("no sub-sampling performed:")
+          bornem0=quantile(Y[1:N/2], probs = c(.2, .8))
+          cat("range m0: ",bornem0,"\n")
+
+          #--- par dichotomie
+          # borne basse
+          N_etape_EM=10
+          m0_tmp=m0_low=bornem0[1]
           m_tmp=m0_tmp
           p_tmp=0.5
-          diff_m0=diff_mm=diff_p=2
+          diff_m0=diff_mm=diff_p=20
           k=1
-          #for (k in 1:N_etape_EM){
-          while (diff_m0 > 0.5 && diff_p > 0.0001){
+          while (diff_m0 > 0.5 && diff_p > 0.0001 && diff_mm > 2){
+            print("ici boucle low")
+            print(k)
             Res_EM=KBO_EM(param=list(mm=m_tmp,
                                      pp=p_tmp,
                                      m0=m0_tmp,
@@ -328,39 +311,151 @@ kfino_fit<-function(datain,Tvar,Yvar,
                           kappaOpt=kappaOpt, Y=Y,Tps=Tps,N=N,dix=dix)
             diff_m0=abs(m0_tmp - Res_EM$m0[[1]])
             diff_p=abs(p_tmp - Res_EM$pp)
+            diff_mm=abs(m_tmp - Res_EM$mm[[1]])
             k<-k+1
-
+            
             m0_tmp=Res_EM$m0[[1]]
             m_tmp=Res_EM$mm[[1]]
             p_tmp=Res_EM$pp
+            print(Res_EM$likelihood)
             if (k==N_etape_EM) break
           }
-          V=Res_EM$likelihood
-
-          if (V > Vopt){
-            Vopt=V
-            m0opt=m0_tmp
-            mmopt=m_tmp
-            popt=p_tmp
+          Vopt_low=Res_EM$likelihood
+          m0opt_low<-Res_EM$m0[[1]]
+          mmopt_low<-Res_EM$mm[[1]]
+          popt_low<-Res_EM$pp
+          
+          # borne haute
+          N_etape_EM=10
+          m0_tmp=m0_up=bornem0[2]
+          m_tmp=m0_tmp
+          p_tmp=0.5
+          diff_m0=diff_mm=diff_p=20
+          k=1
+          while (diff_m0 > 0.5 && diff_p > 0.0001 && diff_mm > 2){
+            print("ici boucle up")
+            print(k)
+            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)
+            diff_m0=abs(m0_tmp - Res_EM$m0[[1]])
+            diff_p=abs(p_tmp - Res_EM$pp)
+            diff_mm=abs(m_tmp - Res_EM$mm[[1]])
+            k<-k+1
+            
+            m0_tmp=Res_EM$m0[[1]]
+            m_tmp=Res_EM$mm[[1]]
+            p_tmp=Res_EM$pp
+            print(Res_EM$likelihood)
+            if (k==N_etape_EM) break
           }
-        }
-
-        print("Optimised parameters with EM method: ")
-        cat("Optimized m0: ",m0opt,"\n")
-        cat("Optimized pp: ",popt,"\n")
-        print("-------:")
-        resultat=KBO_known(param=list(mm=mmopt,
-                                      pp=popt,
-                                      m0=m0opt,
-                                      aa=aa,
-                                      expertMin=expertMin,
-                                      expertMax=expertMax,
-                                      sigma2_m0=sigma2_m0,
-                                      sigma2_mm=sigma2_mm,
-                                      sigma2_pp=sigma2_pp,
-                                      K=K),
-                           threshold=threshold,Y=Y,Tps=Tps,N=N,kappa=kappa)
-
+          Vopt_up=Res_EM$likelihood
+          m0opt_up<-Res_EM$m0[[1]]
+          mmopt_up<-Res_EM$mm[[1]]
+          popt_up<-Res_EM$pp
+          
+          diff_m0range<-abs(m0opt_up - m0opt_low)
+          diff_pprange<-abs(popt_up - popt_low)
+          diff_mmrange<-abs(mmopt_up - mmopt_low)
+          
+          # test quasi equality
+          while(diff_m0range > 0.5 && diff_pprange > 0.0001 && diff_mmrange > 2){
+            m0_tmp=m0_med=(m0_up+m0_low)/2
+            m_tmp=m0_tmp
+            p_tmp=0.5
+            diff_m0=diff_mm=diff_p=20
+            k=1
+            while (diff_m0 > 0.5 && diff_p > 0.0001 && diff_mm > 2){
+               print(k)
+              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)
+              diff_m0=abs(m0_tmp - Res_EM$m0[[1]])
+              diff_p=abs(p_tmp - Res_EM$pp)
+              diff_mm=abs(m_tmp - Res_EM$mm[[1]])
+              k<-k+1
+              
+              m0_tmp=Res_EM$m0[[1]]
+              m_tmp=Res_EM$mm[[1]]
+              p_tmp=Res_EM$pp
+              if (k==N_etape_EM) break
+            }
+            
+            if (Vopt_up < Vopt_low){
+              Vopt_up<-Res_EM$likelihood
+              m0opt_up<-Res_EM$m0[[1]]
+              mmopt_up<-Res_EM$mm[[1]]
+              popt_up<-Res_EM$pp
+              m0_up<-m0_med
+            } else {
+              Vopt_low<-Res_EM$likelihood
+              m0opt_low<-Res_EM$m0[[1]]
+              mmopt_low<-Res_EM$mm[[1]]
+              popt_low<-Res_EM$pp
+              m0_low<-m0_med
+            }
+            
+            diff_m0range<-abs(m0opt_up - m0opt_low)
+            diff_pprange<-abs(popt_up - popt_low)
+            
+          }
+          
+          # challenger
+          m0opt=quantile(Y[1:N/4], probs = c(.5))
+          mmopt=quantile(Y[(3*N/4):N], probs = c(.5))
+          popt=0.5
+          
+          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
+          
+          if (Vopt_low > Vopt){
+            m0opt<-m0opt_low
+            mmopt<-mmopt_low
+            popt<-popt_low
+          }
+          
+          print("Optimised parameters with EM method: ")
+          cat("Optimized m0: ",m0opt,"\n")
+          cat("Optimized mm: ",mmopt,"\n")
+          cat("Optimized pp: ",popt,"\n")
+          print("-------:")
+          resultat=KBO_known(param=list(mm=mmopt,
+                                        pp=popt,
+                                        m0=m0opt,
+                                        aa=aa,
+                                        expertMin=expertMin,
+                                        expertMax=expertMax,
+                                        sigma2_m0=sigma2_m0,
+                                        sigma2_mm=sigma2_mm,
+                                        sigma2_pp=sigma2_pp,
+                                        K=K),
+                             threshold=threshold,Y=Y,Tps=Tps,N=N,kappa=kappa)
+          
       } else if (method == "ML"){
         print("-------:")
         print("Optimisation of initial parameters with ML method - result:")