diff --git a/R/kfino.R b/R/kfino.R index 4821c6e1a596b75b5e650b70cd192169105f688f..ad149a577d3927eae8b31b0bebb0d65485f73e49 100644 --- a/R/kfino.R +++ b/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")