As previously stated, estimation of the hyperparameters \(\hat{r},\hat{\mathbf {\omega }},\hat{\sigma },\hat{\gamma }\)  is conducted by maximizing the marginal likelihood function. Normally, in hierarchical modeling, it is hard to obtain closed-forms for the estimators of the hyperparameters. Thus, numerical optimization is applied to search for optimal hyperparameters. Here, we need to approximate the Hessian matrix at each iteration, which is achieved by applying the Quasi-Newton optimization method to the approximation of the Hessian matrix. The Davidon-Fletcher-Powell (DFP) [1]} formula was proposed as a method to approximate the Hessian matrix, however, this method may not be numerically stable in high-dimensional hyperparameter searches. Rather than solving a linear system of equations as in the DFP the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [2]} was proposed to obtain the search direction using matrix multiplication. The result of this is that the BFGS is an iterative method which is more numerically stable than the DFP, and has very effective “self-correcting properties” [3]} and, thus, performs better for solving unconstrained nonlinear optimization problems such as maximizing the log marginal likelihood presented here.
