R/step3_fit_p.R

In binaryRL: Reinforcement Learning Tools for Two-Alternative Forced Choice Tasks

#' @title 
#' Step 3: Optimizing parameters to fit real data
#' 
#' @description
#' This function is designed to fit the optimal parameters of black-box
#'  functions (models) to real-world data. Provided that the black-box
#'  function adheres to the specified interface
#'  (demo:  
#'     \code{\link[binaryRL]{TD}}, 
#'     \code{\link[binaryRL]{RSTD}}, 
#'     \code{\link[binaryRL]{Utility}}
#'  )
#'  , this function can employ the various optimization algorithms detailed 
#'  below to find the best- fitting parameters for your model.
#'  
#'  The function provides several optimization algorithms:
#'   \itemize{
#'     \item 1. L-BFGS-B (from \code{stats::optim})
#'     \item 2. Simulated Annealing (\code{GenSA::GenSA})
#'     \item 3. Genetic Algorithm (\code{GA::ga})
#'     \item 4. Differential Evolution (\code{DEoptim::DEoptim})
#'     \item 5. Particle Swarm Optimization (\code{pso::psoptim})
#'     \item 6. Bayesian Optimization (\code{mlrMBO::mbo})
#'     \item 7. Covariance Matrix Adapting Evolutionary Strategy (\code{cmaes::cma_es})
#'     \item 8. Nonlinear Optimization (\code{nloptr::nloptr})
#'   }
#' 
#'  For more information, please refer to the homepage of this package:
#'  \url{https://yuki-961004.github.io/binaryRL/}
#'  
#' @note
#' While both \code{fit_p} and \code{rcv_d} utilize the same underlying
#' \code{optimize_para} function to find optimal parameters, they play
#' distinct and sequential roles in the modeling pipeline.
#'
#' The key differences are as follows:
#' \enumerate{
#'   \item{\strong{Purpose and Data Source}: \code{rcv_d} should always be
#'     performed before \code{fit_p}. Its primary role is to validate a
#'     model's stability by fitting it to synthetic data generated by the
#'     model itself. This process, known as parameter recovery, ensures the
#'     model is well-behaved. In contrast, \code{fit_p} is used in the
#'     subsequent stage to fit the validated model to real experimental
#'     data.}
#'
#'   \item{\strong{Estimation Method}: \code{rcv_d} does not include an
#'     \code{estimate} argument. This is because the synthetic data is
#'     generated from known "true" parameters, which are drawn from
#'     pre-defined distributions (typically uniform for most parameters and
#'     exponential for the inverse temperature). Since the ground truth is
#'     known, a hierarchical estimation (MAP) is not applicable. The
#'     \code{fit_p} function, however, requires this argument to handle
#'     real data where the true parameters are unknown.}
#'
#'   \item{\strong{Policy Setting}: In \code{fit_p}, the \code{policy}
#'     setting has different effects: "on-policy" is better for learning
#'     choice patterns, while "off-policy" yields more accurate parameter
#'     estimates. For \code{rcv_d}, the process defaults to an "off-policy"
#'     approach because its main objectives are to verify if the true
#'     parameters can be accurately recovered and to assess whether competing
#'     models are distinguishable, tasks for which off-policy estimation is
#'     more suitable.}
#' }
#'  
#' @param policy [string]
#' 
#' Specifies the learning policy to be used.
#' This determines how the model updates action values based on observed or
#'   simulated choices. It can be either \code{"off"} or \code{"on"}.
#'   
#'  \itemize{
#'   \item {
#'    \strong{Off-Policy (Q-learning): }
#'    This is the most common approach for modeling
#'     reinforcement learning in Two-Alternative Forced Choice (TAFC) tasks.
#'     In this mode, the model's goal is to learn the underlying value of
#'     each option by observing the human participant's behavior. It achieves
#'     this by consistently updating the value of the option that the
#'     human actually chose. The focus is on understanding the value 
#'     representation that likely drove the participant's decisions.
#'   }
#'   \item {
#'    \strong{Off-Policy (SARSA): }
#'    In this mode, the target policy and the behavior policy are identical. 
#'     The model first computes the selection probability for each option based 
#'     on their current values. Critically, it then uses these probabilities to 
#'     sample its own action. The value update is then performed on the action 
#'     that the model itself selected. This approach focuses more on directly 
#'     mimicking the stochastic choice patterns of the agent, rather than just 
#'     learning the underlying values from a fixed sequence of actions.
#'   }
#'  }
#'  
#' default: \code{policy = "off"}
#'  
#' @param estimate [string] 
#' 
#'   Estimation method. Can be either \code{"MLE"} or \code{"MAP"}.
#'   \itemize{
#'     \item{\code{"MLE"}: (Default) Maximum Likelihood Estimation. This
#'       method finds the parameter values that maximize the log-likelihood
#'       of the data. A higher log-likelihood indicates that the parameters
#'       provide a better explanation for the observed human behavior. In
#'       other words, data simulated using these parameters would most
#'       closely resemble the actual human data. This method does not
#'       consider any prior information about the parameters.}
#'       
#'     \item{\code{"MAP"}: Maximum A Posteriori Estimation. This method
#'       finds the parameter values that maximize the posterior probability.
#'       It is an iterative process based on the Expectation-Maximization
#'       (EM) framework.
#'       
#'       \itemize{
#'         \item{\strong{Initialization}: The process begins by assuming a
#'           uniform distribution as the prior for each parameter, making the
#'           initial log-prior zero. The first optimization is thus
#'           equivalent to MLE.}
#'         \item{\strong{Iteration}: After finding the best parameters for
#'           all subjects, the algorithm assesses the actual distribution of
#'           each parameter and fits a normal distribution to it. This
#'           fitted distribution becomes the new empirical prior.}
#'         \item{\strong{Re-estimation}: The parameters are then re-optimized
#'           to maximize the updated posterior probability.}
#'         \item{\strong{Convergence}: This cycle repeats until the posterior
#'           probability converges or the maximum number of iterations
#'           (specified by \code{iteration_g}) is reached.}
#'       }
#'       Using this method requires that the \code{priors} argument
#'       be specified to define the initial prior distributions.
#'     }
#'    }
#' 
#' default: \code{estimate = "MLE"}
#'  
#' @param data [data.frame] 
#' 
#' This data should include the following mandatory columns: 
#'  \itemize{
#'    \item \code{sub} "Subject"
#'    \item \code{time_line} "Block" "Trial"
#'    \item \code{L_choice} "L_choice"
#'    \item \code{R_choice} "R_choice"
#'    \item \code{L_reward} "L_reward"
#'    \item \code{R_reward} "R_reward"
#'    \item \code{sub_choose} "Sub_Choose"
#'  }
#'  
#' @param id [CharacterVector]
#' 
#'  A vector specifying the subject ID(s) for which parameters should be
#'   fitted. The function will process only the subjects provided in this
#'   vector.
#'
#'  To fit all subjects, you can either explicitly set the argument as
#'   \code{id = unique(data$Subject)} or leave it as the default
#'   (\code{id = NULL}). Both approaches will direct the function to fit
#'   parameters for every unique subject in the dataset.
#'
#'  It is strongly recommended to avoid using simple numeric sequences like
#'   \code{id = 1:4}. This practice can lead to errors if subject IDs are
#'   stored as strings (e.g., subject four is stored as "004") or are not
#'   sequentially numbered.
#'  
#'  default: \code{id = NULL}
#'  
#' @param n_trials [integer]
#' 
#' Represents the total number of trials a single subject experienced
#'  in the experiment. If this parameter is kept at its default value
#'  of \code{NULL}, the program will automatically detect how many trials
#'  a subject experienced from the provided data. This information
#'  is primarily used for calculating model fit statistics such as
#'  AIC (Akaike Information Criterion) and BIC (Bayesian Information
#'  Criterion).
#'  
#'  default: \code{n_trials = NULL}
#' 
#' @param funcs [CharacterVector]
#' 
#' A character vector containing the names of all user-defined functions
#'  required for the computation. When parallel computation is enabled
#'  (i.e., \code{nc > 1}), user-defined models and their custom functions 
#'  might not be automatically accessible within the parallel environment.
#'
#' Therefore, if you have created your own reinforcement learning model
#'  that modifies the package's default six default functions 
#'  (default functions: 
#'     \code{util_func = \link[binaryRL]{func_gamma}}, 
#'     \code{rate_func = \link[binaryRL]{func_eta}}, 
#'     \code{expl_func = \link[binaryRL]{func_epsilon}},
#'     \code{bias_func = \link[binaryRL]{func_pi}},
#'     \code{prob_func = \link[binaryRL]{func_tau}},
#'     \code{loss_func = \link[binaryRL]{func_logl}}
#'  ), 
#'  you must explicitly provide the names of your custom functions as a 
#'  vector here.
#' 
#' @param model_name [List] 
#' 
#' The name of fit modals
#' 
#' e.g. \code{model_name = c("TD", "RSTD", "Utility")}
#' 
#' @param fit_model [List] 
#' 
#' A collection of functions applied to fit models to the data.
#' 
#' e.g. \code{fit_model = list(binaryRL::TD, binaryRL::RSTD, binaryRL::Utility)}
#' 
#' @param lower [List] 
#' 
#' The lower bounds for model fit models
#' 
#' e.g. \code{lower = list(c(0, 0), c(0, 0, 0), c(0, 0, 0))}
#' 
#' @param upper [List] 
#' 
#' The upper bounds for model fit models
#' 
#' e.g. \code{upper = list(c(1, 10), c(1, 1, 10), c(1, 1, 10))}
#' 
#' @param priors [List]
#' 
#'  A list specifying the prior distributions for the model parameters.
#'   This argument is mandatory when using \code{estimate = "MAP"}.
#'   There are two primary scenarios for its use:
#'
#'   \strong{1. Static MAP Estimation (Non-Hierarchical)}
#'   This approach is used when you have a strong, pre-defined belief about
#'   the parameter priors and do not want the model to update them iteratively.
#'   \describe{
#'     \item{Configuration:}{
#'       \itemize{
#'         \item{Set \code{estimate = "MAP"}.}
#'         \item{Provide a \code{list} defining your confident prior
#'           distributions.}
#'         \item{Keep \code{iteration_g = 0} (the default).}
#'       }
#'     }
#'     \item{Behavior:}{ The algorithm maximizes the posterior
#'       probability based solely on your specified priors. It will not
#'       use the EM (Expectation-Maximization) framework to learn new
#'       priors from the data.}
#'   }
#'
#'   \strong{2. Hierarchical Bayesian Estimation via EM}
#'   This approach is used to let the model learn the group-level (hierarchical)
#'   prior distributions directly from the data.
#'   \describe{
#'     \item{Configuration:}{
#'       \itemize{
#'         \item{Set \code{estimate = "MAP"}.}
#'         \item{Specify a weak or non-informative initial prior, such as a
#'           uniform distribution for all parameters.}
#'         \item{Set \code{iteration_g} to a value greater than 0.}
#'       }
#'     }
#'     \item{Behavior:}{ With a uniform prior, the initial
#'       log-posterior equals the log-likelihood, making the first estimation
#'       step equivalent to MLE. The algorithm then initiates the EM procedure:
#'       it iteratively assesses the actual parameter distribution across all
#'       subjects and updates the group-level priors. This cycle continues
#'       until the posterior converges or \code{iteration_g} is reached.}
#'   }
#' 
#'  default: \code{priors = NULL}
#' 
#' @param tolerance [double] 
#' 
#' Convergence threshold for MAP estimation. If the change in
#'  log posterior probability between iterations is smaller than this value, the
#'  algorithm is considered to have converged and the program will stop.
#' 
#' default: \code{tolerance = 0.001}
#' 
#' @param initial_params [NumericVector]
#' 
#' Initial values for the free parameters that the optimization algorithm will
#'  search from. These are primarily relevant when using algorithms that require
#'  an explicit starting point, such as \code{L-BFGS-B}. If not specified,
#'  the function will automatically generate initial values close to zero.
#'  
#'  default: \code{initial_params = NA}.
#'
#' @param initial_size [integer]
#' 
#' This parameter corresponds to the population size in genetic 
#'  algorithms (\code{GA}). It specifies the number of initial candidate
#'  solutions that the algorithm starts with for its evolutionary search.
#'  This parameter is only required for optimization algorithms that operate on
#'  a population, such as \code{GA} or \code{DEoptim}. 
#'  
#'  default: \code{initial_size = 50}.
#'  
#' @param iteration_i [integer] 
#' 
#' The number of iterations the optimization algorithm will perform
#'  when searching for the best-fitting parameters during the fitting
#'  phase. A higher number of iterations may increase the likelihood of 
#'  finding a global optimum but also increases computation time.
#'  
#'  default: \code{iteration_i = 10}.
#'  
#' @param iteration_g [integer]
#' 
#' The maximum number of iterations for the Expectation-Maximization (EM)
#'  based MAP estimation. The algorithm will stop once this iteration
#'  count is reached, even if the change in the log-posterior value has
#'  not yet fallen below the \code{tolerance} threshold.
#' 
#'  default: \code{iteration_g = 0}.
#' 
#' @param seed [integer] 
#' 
#' Random seed. This ensures that the results are 
#'  reproducible and remain the same each time the function is run. 
#'  
#' default: \code{seed = 123}
#'  
#' @param nc [integer]
#' 
#' Number of cores to use for parallel processing. Since fitting
#' optimal parameters for each subject is an independent task,
#' parallel computation can significantly speed up the fitting process:
#' \itemize{
#'   \item \strong{nc = 1}: The fitting proceeds sequentially.
#'   Parameters for one subject are fitted completely before moving
#'   to the next subject.
#'   \item \strong{nc > 1}: The fitting is performed in parallel
#'   across subjects. For example, if \code{nc = 4}, the algorithm will
#'   simultaneously fit data for four subjects. Once these are complete,
#'   it will proceed to fit the next batch of subjects (e.g., subjects
#'   5-8), and so on, until all subjects are processed.
#' }
#' 
#'  \code{default: nc = 1}
#'  
#' @param algorithm [string] 
#' 
#' Choose an algorithm package from
#'  \code{L-BFGS-B}, \code{GenSA},\code{GA},\code{DEoptim},\code{PSO},
#'  \code{Bayesian}, \code{CMA-ES}.
#'  
#' In addition, any algorithm from the \code{nloptr} package is also
#'  supported. If your chosen \code{nloptr} algorithm requires a local search,
#'  you need to input a character vector. The first element represents
#'  the algorithm used for global search, and the second element represents
#'  the algorithm used for local search.
#'
#' @returns 
#' The optimal parameters found by the algorithm for each subject,
#'  along with the model fit calculated using these parameters.
#'  This is returned as an object of class \code{binaryRL} containing results
#'  for all subjects with all models.
#'  
#' @examples
#' \dontrun{
#' comparison <- binaryRL::fit_p(
#'   data = binaryRL::Mason_2024_G2,
#'   id = unique(binaryRL::Mason_2024_G2$Subject),
#' #+-----------------------------------------------------------------------------+#
#' #|----------------------------- black-box function ----------------------------|#
#'   #funcs = c("your_funcs"),
#'   policy = c("off", "on"),
#'   fit_model = list(binaryRL::TD, binaryRL::RSTD, binaryRL::Utility),
#'   model_name = c("TD", "RSTD", "Utility"),
#' #|--------------------------------- estimate ----------------------------------|#
#'   estimate = c("MLE", "MAP"),
#' #|------------------------------------ MLE ------------------------------------|#
#'   lower = list(c(0, 0), c(0, 0, 0), c(0, 0, 0)),
#'   upper = list(c(1, 10), c(1, 1, 10), c(1, 1, 10)),
#' #|------------------------------------ MAP ------------------------------------|#
#'   priors = list(
#'     list(
#'       eta = function(x) {stats::dunif(x, min = 0, max = 1, log = TRUE)}, 
#'       tau = function(x) {stats::dexp(x, rate = 1, log = TRUE)}
#'     ), 
#'     list(
#'       eta = function(x) {stats::dunif(x, min = 0, max = 1, log = TRUE)}, 
#'       eta = function(x) {stats::dunif(x, min = 0, max = 1, log = TRUE)}, 
#'       tau = function(x) {stats::dexp(x, rate = 1, log = TRUE)}
#'     ), 
#'     list(
#'       eta = function(x) {stats::dunif(x, min = 0, max = 1, log = TRUE)}, 
#'       gamma = function(x) {stats::dunif(x, min = 0, max = 1, log = TRUE)}, 
#'       tau = function(x) {stats::dexp(x, rate = 1, log = TRUE)}
#'     )
#'   ),
#' #|----------------------------- iteration number ------------------------------|#
#'   iteration_i = 10,
#'   iteration_g = 10,
#' #|-------------------------------- algorithms ---------------------------------|#
#'   nc = 1,                 # <nc > 1>: parallel computation across subjects
#'   # Base R Optimization
#'   algorithm = "L-BFGS-B"  # Gradient-Based (stats)
#' #|-----------------------------------------------------------------------------|#
#'   # Specialized External Optimization
#'   #algorithm = "GenSA"    # Simulated Annealing (GenSA)
#'   #algorithm = "GA"       # Genetic Algorithm (GA)
#'   #algorithm = "DEoptim"  # Differential Evolution (DEoptim)
#'   #algorithm = "PSO"      # Particle Swarm Optimization (pso)
#'   #algorithm = "Bayesian" # Bayesian Optimization (mlrMBO)
#'   #algorithm = "CMA-ES"   # Covariance Matrix Adapting (cmaes)
#' #|-----------------------------------------------------------------------------|#
#'   # Optimization Library (nloptr)
#'   #algorithm = c("NLOPT_GN_MLSL", "NLOPT_LN_BOBYQA")
#' #|-------------------------------- algorithms ---------------------------------|#
#' #################################################################################
#' )
#'
#' result <- dplyr::bind_rows(comparison)
#'
#' # Ensure the output directory exists before writing
#' if (!dir.exists("../OUTPUT")) {
#'   dir.create("../OUTPUT", recursive = TRUE)
#' }
#' 
#' write.csv(result, "../OUTPUT/result_comparison.csv", row.names = FALSE)
#' }
#' 
fit_p <- function(
  policy = "off",
  estimate = "MLE",  
  
  data,
  id = NULL,
  n_trials = NULL,
  
  funcs = NULL,
  model_name = c("TD", "RSTD", "Utility"),
  fit_model = list(binaryRL::TD, binaryRL::RSTD, binaryRL::Utility),
  lower = list(c(0, 0), c(0, 0, 0), c(0, 0, 0)),
  upper = list(c(1, 1), c(1, 1, 1), c(1, 1, 1)),
  
  priors = NULL,
  tolerance = 0.001,
  
  iteration_i = 10,
  iteration_g = 0,
  
  initial_params = NA,
  initial_size = 50,
  seed = 123,
  
  nc = 1,
  algorithm
){
################################# [ info ] #####################################   
  
  # 事前准备. 探测信息
  info <- suppressWarnings(suppressMessages(detect_information(data = data)))
  
  if (is.null(n_trials)) {
    n_trials <- info[["n_trials"]]
  }
  
  if (is.null(id)) {
    id <- info[["all_ids"]]
  }
  
  # 创建空list用于储存结果
  model_comparison <- list()
  model_result <- list()

################################ [ threads ] ###################################  
  
  sys <- Sys.info()[["sysname"]]
  
  if (nc == 1) {
    future::plan(future::sequential)
  } 
  # Windows
  else if (sys == "Windows") {
    future::plan(future::multisession, workers = nc)
  } 
  # macOS
  else if (sys == "Darwin") {  
    future::plan(future::multisession, workers = nc)
  } 
  # Linux
  else if (sys == "Linux") {
    future::plan(future::multicore, workers = nc)
  }
  
  doFuture::registerDoFuture()

############################# [ for loop models] ###############################
    
  # 每个model都使用初始的priors
  params_priors <- list()
  
  for (i in 1:length(fit_model)){
    
    # 编译对象函数
    obj_func <- compiler::cmpfun(fit_model[[i]])
    
    message(paste0(
      "\n", 
      "Fitting Model: ", model_name[i], 
      "\n"
    ))
    
    n_subjects <- length(id)
    n_params <- length(lower[[i]])

########################### [ foreach loop subs ] ##############################    
    
    # 进度条
    progressr::handlers(progressr::handler_txtprogressbar)
    
    progressr::with_progress({
      
      # 进度条与此时运行的被试数绑定
      p <- progressr::progressor(steps = n_subjects)
      
      # doRNG保证种子不变
      doRNG::registerDoRNG(seed = seed)
      
      # 初始化(定义)foreach中的j
      j <- NA

################################# [ MLE ] ######################################
      
      # 抑制每个线程加载包时的信息
      suppressMessages({
        model_result <- foreach::foreach(
          j = 1:n_subjects, .combine = rbind,
          .packages = c("binaryRL"),
          .export = funcs
        ) %dorng% {

          binaryRL_res <- binaryRL::optimize_para(
            policy = policy,
            estimate = estimate,
            
            data = data,
            id = id[j],
            n_trials = n_trials,
            n_params = n_params,
            
            obj_func = obj_func,
            lower = lower[[i]],
            upper = upper[[i]],
            priors = priors[[i]],
            
            initial_params = initial_params,
            initial_size = initial_size,
            
            iteration = iteration_i,
            seed = seed,
            algorithm = algorithm
          )
          
          result_j <- data.frame(
            fit_model = model_name[i],
            Subject = id[j],
            ACC = binaryRL_res$acc,
            LogL = -binaryRL_res$ll,
            LogPr = -binaryRL_res$lpr,
            LogPo = -binaryRL_res$lpo,
            AIC = binaryRL_res$aic,
            BIC = binaryRL_res$bic
          )
          
          for (k in 1:n_params) {
            result_j[1, k + 8] <- binaryRL_res$output[k]
            names(result_j)[k + 8] <- paste0("param_", k)
          }
          
          # 在foreach循环内更新进度条
          p() 
          return(result_j)
        }
      })
    })

################################# [ MAP ] ######################################
        
    # Expectation-Maximization Algorithm
    if (estimate == "MAP") {
      
      # 基于上面第一次的结果更新先验概率
      params_priors[[i]] <- update_priors(
        model_result = model_result,
        priors = priors[[i]], 
        n_params = n_params,
        param_prefix = "param_"
      )
      
      # 计算第一次拟合时的似然后验值
      LogPo <- sum(model_result$LogPo)
      
      # 因为之后使用while循环
      # 如果iteration_g > 0, 则说明用户希望进行层级贝叶斯, 至少进行一次EM
      if (iteration_g > 0) {
        delta_LogPo <- tolerance + 1
        message(
          paste0(
            "Starting Expectation-Maximization Algorithm", "\n",
            "Log-Posterior Probability: ", round(LogPo, 2)
          )
        )
      }
      # 如果iteration_g == 0, 则意味着被试想只进行MAP或MLE, 而不是层级贝叶斯
      else {
        delta_LogPo <- 0
      }
      
      # 初始化后验迭代次数
      iteration <- 0
      
      while (abs(delta_LogPo) > tolerance) {
        
        # 进度条
        progressr::handlers(progressr::handler_txtprogressbar)
        
        progressr::with_progress({
          
          # 进度条与此时运行的被试数绑定
          p <- progressr::progressor(steps = n_subjects)
          
          # doRNG保证种子不变
          doRNG::registerDoRNG(seed = seed)
          
          # 初始化(定义)foreach中的j (无意义, 仅为通过R CMD check)
          j <- NA
          
          # 抑制每个线程加载包时的信息
          suppressMessages({
            model_result <- foreach::foreach(
              j = 1:n_subjects, .combine = rbind,
              .packages = c("binaryRL"),
              .export = funcs
            ) %dorng% {
              
              binaryRL_res <- binaryRL::optimize_para(
                policy = policy,
                estimate = estimate,
                
                data = data,
                id = id[j],
                n_trials = n_trials,
                n_params = n_params,
                
                obj_func = obj_func,
                lower = lower[[i]],
                upper = upper[[i]],
                priors = params_priors[[i]],
                
                initial_params = initial_params,
                initial_size = initial_size,
                
                iteration = iteration_i,
                seed = seed,
                algorithm = algorithm
              )
              
              result_j <- data.frame(
                fit_model = model_name[i],
                Subject = id[j],
                ACC = binaryRL_res$acc,
                LogL = -binaryRL_res$ll,
                LogPr = -binaryRL_res$lpr,
                LogPo = -binaryRL_res$lpo,
                AIC = binaryRL_res$aic,
                BIC = binaryRL_res$bic
              )
              
              for (k in 1:n_params) {
                result_j[1, k + 8] <- binaryRL_res$output[k]
                names(result_j)[k + 8] <- paste0("param_", k)
              }
              
              # 在foreach循环内更新进度条
              p() 
              return(result_j)
            }
          })
        })
        
        # 计算此时的似然后验变化量
        delta_LogPo <- sum(model_result$LogPo) - LogPo
        
        # 更新似然后验
        LogPo <- sum(model_result$LogPo)
        
        message(
          paste0(
            "Log-Posterior Probability: ", round(LogPo, 2),
            ", ",
            "\u0394: ", sign_numbers(delta_LogPo), round(delta_LogPo, 2)
          )
        )
        
        # 修改先验分布
        params_priors[[i]] <- update_priors(
          model_result = model_result,
          priors = params_priors[[i]], 
          n_params = n_params,
          param_prefix = "param_"
        )
        
        # EM也不会无限制跑下去
        iteration <- iteration + 1
        if (iteration >= iteration_g) {
          break 
        }
      }
    }
    
    # 每次将一个模型的结果输出, 然后跑下一个模型的拟合
    model_comparison[[i]] <- model_result
  }
  
  # 取消注册的线程
  future::plan(future::sequential)
  
  result <- model_comparison
  
  return(result)
}