update KBO_EM - while loop

R/kfino.R

@@ -7,8 +7,8 @@
 #' @param doOptim logical, if TRUE optimisation of the initial parameters
 #'                default TRUE
 #' @param method character, the method used to optimize the initial parameters:
-#'               Expectation-Maximization algorithm `"EM"` or Maximization 
-#'               Likelihhod `"ML"`, default `"ML"` 
+#'               Expectation-Maximization algorithm `"EM"` or Maximization
+#'               Likelihhod `"ML"`, default `"ML"`
 #' @param threshold numeric, threshold to qualify an observation as outlier
 #'        according to the label_pred, default 0.5
 #' @param kappa numeric, truncation setting, default 10
@@ -40,7 +40,7 @@
 #' range) using quantile of the Y distribution (varying between 0.2 and 0.8 for
 #' m0 and 0.5 for mm). pp is a sequence varying between 0.5 and 0.7. A
 #' sub-sampling is performed to speed the algorithm if the number of possible
-#' observations studied is greater than 500. Optimization is performed using 
+#' observations studied is greater than 500. Optimization is performed using
 #' `"EM"` or `"ML"` methods.
 #'
 #' @importFrom stats dnorm quantile na.omit
@@ -79,14 +79,14 @@
 #'               Tvar="dateNum",Yvar="Poids",
 #'               doOptim=TRUE,method="ML",param=param1)
 #' Sys.time() - t0
-#' 
+#'
 #' # --- With Optimisation on initial parameters - EM method
 #' t0 <- Sys.time()
 #' resu1b<-kfino_fit(datain=spring1,
 #'               Tvar="dateNum",Yvar="Poids",
 #'               doOptim=TRUE,method="EM",param=param1)
 #' Sys.time() - t0
-#' 
+#'
 #' # --- Without Optimisation on initial parameters
 #' t0 <- Sys.time()
 #' param2<-list(m0=41,
@@ -290,7 +290,7 @@ kfino_fit<-function(datain,Tvar,Yvar,
         cat("range m0: ",bornem0,"\n")
         cat("initial m0opt: ",m0opt,"\n")
         popt=0.5
-        
+
         # 1er tour
         Vopt=KBO_EM(param=list(m0=m0opt,
                                 mm=mmopt,
@@ -303,7 +303,7 @@ kfino_fit<-function(datain,Tvar,Yvar,
                                 sigma2_pp=sigma2_pp,
                                 K=K),
                      Y=Y,Tps=Tps,N=N,dix=dix,kappaOpt=kappaOpt)$likelihood
-        
+
         # boucle
         N_etape_EM=10
         t=1
@@ -312,7 +312,11 @@ kfino_fit<-function(datain,Tvar,Yvar,
           t=t+1
           m_tmp=m0_tmp
           p_tmp=0.5
-          for (k in 1:N_etape_EM){
+          diff_m0=diff_mm=diff_p=2
+          k=1
+          #for (k in 1:N_etape_EM){
+          while (k <= N_etape_EM || (diff_m0 > 0.5 && diff_p > 0.0001) ){
+            print("passage")
             Res_EM=KBO_EM(param=list(mm=m_tmp,
                                      pp=p_tmp,
                                      m0=m0_tmp,
@@ -324,22 +328,28 @@ kfino_fit<-function(datain,Tvar,Yvar,
                                      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)
+            print(diff_m0)
+            print(diff_p)
+            k<-k+1
+
             m0_tmp=Res_EM$m0[[1]]
             m_tmp=Res_EM$mm[[1]]
             p_tmp=Res_EM$pp
           }
           V=Res_EM$likelihood
-          
+
           print(m0_tmp)
           if (V > Vopt){
-            print("EM converged.")
+            #print("EM converged.")
             Vopt=V
             m0opt=m0_tmp
             mmopt=m_tmp
             popt=p_tmp
           }
         }
-        
+
         print("Optimised parameters with EM method: ")
         cat("Optimized m0: ",m0opt,"\n")
         cat("Optimized pp: ",popt,"\n")
@@ -355,7 +365,7 @@ kfino_fit<-function(datain,Tvar,Yvar,
                                       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:")
@@ -363,13 +373,13 @@ 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
-        
+
         Vopt=KBO_L(list(m0=m0opt,
                         mm=mmopt,
                         pp=popt,
@@ -395,7 +405,7 @@ kfino_fit<-function(datain,Tvar,Yvar,
                            sigma2_pp=sigma2_pp,
                            K=K),
                       Y=Y,Tps=Tps,N=N,dix=dix,kappaOpt=kappaOpt)
-              
+
               if (V > Vopt){
                 Vopt=V
                 m0opt=m0
@@ -405,7 +415,7 @@ kfino_fit<-function(datain,Tvar,Yvar,
             }
           }
         }
-        
+
         print("Optimised parameters: ")
         cat("Optimized m0: ",m0opt,"\n")
         cat("Optimized mm: ",mmopt,"\n")