Isomerization of α-pinene is to produce other monoterpene hydrocarbons that have wide application in the cosmetic, perfume, aromatic polymer, and other domestic chemicals industries  . Thermal isomerization of α-pinene in liquid-phase  , gas-phase  , and supercritical solvents-fluids  is studied to obtain higher production yields and selectivity. Reaction kinetic modeling is a vital tool for the simulation of the thermal isomerization of α-pinene and for industrial reactor design. Reaction modeling consists of two main steps: first is to model the reaction pathways known as a kinetic model that needs an insight understanding of mechanisms in the reaction background and second is to calibrating the kinetic model by using experimental data  .
Kinetic models often involve a number of rate constants to be estimated. The value of these parameters is heavily dependent on the reaction conditions such as temperature, pressure, reaction media, the assumed reaction pathway, and the quantity and quality of the measured experimental data. These are called uncertainty sources and can significantly affect the estimated value of the reaction rates.
Least-square regressions are extensively used to find the best-fitted curve to the experimental data with a given kinetic model. Classical gradient-based optimization methods such as gradient descent and evolutionary algorithms (EA) and Genetic algorithm are applied to minimize the square error between the measured and modeled data points. In an ill-conditioned system of equations, the gradient-based algorithms can trap in local minimums and are usually failed to find the global minimum error  . However, the EA methods result in a global optimum. Kinetic model of the thermal isomerization of α-pinene is studied by Tjoa et al.  who applied local estimator methods for the rate constant estimation and found that these local methods need good initial guesses to be able to find the optimum rate constants. Yermakova et al.  applied gradient- based Gauss-Marquardt iteration method for the kinetic model identification of the thermal isomerization of α-pinene and in order to overcome to the ill-con- ditioning conditions, they applied regularization techniques that can result in a biased estimation of reaction constants. Rodriguez-Fernandez et al.  studied the same reaction system and successfully applied a global optimization based on the Scatter Search method and obtained the point estimate of rate constants and the global minimum of the least-square objective function. However, these approaches result in a point estimate of reaction rates without taking into account the uncertainty. Therefore, other statistical techniques such as analysis of variance (ANOVA), t-distribution tests, Fisher Information Matrix, jackknifing, and bootstrapping methods need to be used as complementary methods to obtain the parameter covariance matrix.
Bayesian inference is a logical statistical framework that applies experimental measured data as evidence points to increase the degree of confidence for the targeted parameters with a given prior knowledge. In the kinetic model calibration, the Bayesian inference can be used to quantify the uncertainty associated with the rate constant values as well as model structural error. Markov chain Monte Carlo (MCMC) methods such as Metropolis-Hastings algorithm are usually used to construct the posterior distribution of the Bayesian inference. This method results in an evidence base confidence interval of the rate constants while the correlation between parameters is taken into consideration. MCMC methods need to take large enough samples from the posterior distribution to result in a converged histogram of parameters. Therefore, applying MCMC methods are computationally expensive. Galagali et al.  proposed an adaptive MCMC to alleviate a portion of this computational burdensome. High-perfor- mance computing can largely enhance the performance of MCMC sampling methods to deal with uncertainties. For example, Alikhani et al.   showed an application of adaptive ODE solution algorithm in parallel-based GPU to accelerate the solution of mass balance rate equations in the kinetic models. They showed that applying these improvements can reduce the overall computation time by up to 50%. Altogether, the nowadays high-speed computers plus high- performance computing algorithms made it more feasible to use MCMC methods in the uncertainty evaluation.
Moles et al.  studied a non-linear biochemical model and applied different optimization method to estimate 36 model parameters and concluded that only stochastic algorithms were able to successfully solve the problem. Sathyamoorthy et al.  applied the Generalized Uncertainty Estimation (GLUE) technique for the uncertainty evaluation of a biological nitrification model. Sun et al.  applied the Bayesian inference with MCMC to obtain rate constants of Phytate (an organic phosphorus compound in agricultural soils) degradation in three different pathways. The Bayesian inference is applied in other studies in different fields for uncertainty analysis and parameter estimation      .
This study is aiming at applying the Bayesian inference to numerically evaluate the uncertainty associated with the reaction rate constants and their correlation matrix for the thermal isomerization of α-pinene over the experimental data obtained from Fuguitt et al.  . The obtained results are compared with the results in other studies obtained by using other methods to demonstrate the strength of proposed method.
2. Materials and Methods
Experimental data of liquid-phase thermal isomerization of α-pinene in 189.5°C obtained by Fuguitt et al.  is presented in Table 1. This reaction produces four products that are presented in weight fraction in eight intervals each with duplicated results. Hunter et al.  proposed a kinetic model for this reaction as it is shown in Figure 1 followed by Equations (1) to (5). In this scheme α-pinene (C1) converts to dipentene (C2) and allo-ocimen (C3) in two parallel paths; dipentene is the end-point component but allo-ocimen yields α- and β-pyronene (C4) and involves in an equilibrium reaction with other dimer compounds (C5).
In this kinetic model, each reaction pathway follows by a first order reaction.
Figure 1. Kinetic model of the thermal isomerization of α-pinene where C1 is α-pinene, C2 is dipenetene, C3 is allo-ocimene, C4 is pyrone, and C5 denotes for other dimer compounds  .
Table 1. Experimental data of thermal isomerization of α-pinene at 189.5 °C  .
The kinetic model can be simplified to show the system of ordinary differential equations (ODEs):
where j denotes each reaction pathway and is the stoichiometery of component i in path j. By introducing as normalized concentration where
is the initial concentration of C1 (α-pinene), the Equation (6) can be rewritten as  :
With the following initial conditions:, and other equal to zero.
Model calibration is determining the numerical value of rate constants () considering the experimental data.
3. Model Calibration
The idea of applying the Bayesian inference in the rate constant evaluation of a kinetic model is based on recent work by Alikhani et al.  , which they proposed a systematic way to apply the Bayesian inference for parameter estimation and uncertainty quantification of a large-sized kinetic model.
Alikhani et al.  introduced a general form of the kinetic model in a vector form:
where f represents the interactions between species by applying mass balance equation while is the concentrations of all the species, is the external forcing input, and is the model parameters. In the batch reactor system there is no external forcing input and therefore the concentration of each species becomes only a function of model parameters and the initial concentrations.
In the normalized form followed by Equation (7), the model output concentrations become only a function of model parameters:
By using this notation, model calibration defines as estimating the values of values approaching a maximum likelihood of modeled () versus experimental concentrations. Assuming a normal distribution for likelihood function with constant error standard deviation (ESD), , the likelihood probability distribution can be define as  :
where and are the individual points of and vectors, respectively, and n is the total number of experimental data. According to Equation (10), is only a function of in a given initial conditions and therefore is equivalent term for the likelihood probability function.
By this assumption (normal distribution for errors and a constant independent standard error), the square error minimization becomes a specific case of maximum likelihood method as minimizing the numerator of the exponential in Equation (11) maximizes the. This is referred to as deterministic approach when the ESD assumes to be constant; otherwise, it is called stochastic approach when ESD treats as an unknown variable resulting in a distribution of valid parameters instead of a point estimate. Therefore, in the stochastic approach the likelihood function can be defined as where is the union of rate constants and the ESD as. Estimating the numerical value of the ESD is referred to as uncertainty quantification.
Stochastic model calibration is to find all set of rate constants that meet the experimental data. In order to find the, we can use the Bayesian inference    :
where is the prior distribution of parameters, is the posterior distribution of parameters that is deemed to be estimated, and is the probability distribution of experimental data that is unknown and makes direct application of the Bayesian inference impossible.
The prior distribution can take a uniform distribution within the lower and upper range of parameters, a normal distribution with its prior mean and standard deviation, or another type of distributions. For example, Alikhani et al.  assumed that log-normal distribution is a better representation of the prior distribution of the parameters in their kinetic model. In this study, to be consistent with the likelihood distribution, all the prior distributions are assumed to be a normal distribution.
MCMC is a class of techniques to sample from multidimensional joint probability distributions     . Random walk Metropolis-Hastings (MH) algorithm   is a subset of MCMC algorithms that is used in this study to estimate the probability histogram of parameters, defines as:
if else (13)
where is a uniform random number between 0 and 1, is the
where is the proposal distribution that randomly proposes by given. If the Metropolis-Hastings uses as random walk Monte Carlo then g must be symmetric  and the acceptance rate simplifies as:
where is the posterior probability of parameter that is unknown. However, by applying Bayesian inference the acceptance rate can be obtained as:
where by given and all the terms in the right-hand side can be calculated. If g (the proposal distribution) is chosen to be a normal distribution then:
where is a random number from the standard normal distribution, is the prior standard deviation of parameters, and is an adjustable factor that can take any positive value less than unity to adjust the convergence speed of the random walk. Convergence in the MCMC is defined as conditions that the change in the standard deviation of the posterior distribution becomes less than threshold when random walk sampling continued.
4. Results and Discussion
To obtain the posterior distribution of the parameters, Equation (13) was implemented in MATLAB and the ODE equation set of Equation (7) was solved by ode45 built-in function. Alikhani et al.  suggested throwing away the initial portion of the random walk samples due to burn-in period. Therefore, after 5000 samples as burn-in samples, the MCMC convergence was examined every 10,000 samples, and it is found that after 180,000 samples the change in the standard deviation of posterior distributions is negligible. Histograms of the five rate constants and the ESDn are shown in Figure 2. The 95% confidence interval of the rate constants is extracted from the posterior distribution and is shown in Table 2.
The deterministic point estimate of the parameters is also obtained by running Genetic algorithm (GA). Comparing the point estimate results shows a very good agreement between the point estimates and the median of the 95% posterior interval of the 5 rate constants. However, the value obtained for the ESD in the deterministic approach with the value of 0.77 is significantly different than its confidence interval of 1.73 - 2.35 obtained by using Bayesian inference application. This is the main difference between the deterministic and the stochastic approaches which the ESD―the representative of the uncertainty―is missing in the deterministic calculations. Therefore, the ESD of the model is not true in the point estimate method.
The ESD of the model is representing both the uncertainty associated with the measurements error and the uncertainty associated with the model structure. It should be noted that for each species in the kinetic model a separate ESD could be considered resulting in a more insight uncertainty assessment. However, adding more parameters to the random walk space needs more experimental data and larger sample size (resulting in higher computational cost) to be converged.
Figure 2. Histograms representing the posterior probability density of the thermal isomerization of α-pinene rate constants (105 min−1) and the error standard deviation.
Rodriguez-Fernandez et al.  studied the same kinetic model for the isomerization of α-pinene and applied a global optimization based on the Scatter
Table 2. The confidence interval of the thermal isomerization of α-pinene rate constants (105 min−1) and the error standard deviation of the kinetic model.
Search method for model calibration. They obtained the point estimates of p1 = 5.9259, p2 = 2.9634, p3 = 2.0473, p4 = 27.449, and p5 = 3.9980 (values are times to 105 to be consistent with the units in Table 2). The results of their study are in a close agreement with the results obtained in this study. However, this was not a surprising point because both studies are applied an evolutionary optimization technique on the same set of experimental data. But what distinguishes two studies is the significance evaluation of the parameter values. Rodriguez-Fernandez et al.  applied an approximation method based on a local linearization of the model output to evaluate the parameter covariance matrix based on the Fisher information matrix. They obtained the value of 0.07195, 0.06519, 0.33328, 2.7657, and 0.9757 for the standard deviation of p1 to p5, respectively; where comparing with the standard deviation of parameters in Table 2, it shows roughly 86%, 76%, 33%, 38%, and 28% higher values. This point is also mentioned by Rodriguez-Fernandez et al.  that the confidence intervals obtained from the Fisher information matrix show a wider confidence interval. It can be concluded that the confidence interval of rate constants obtained by applying Bayesian inference concept is a better representation of uncertainty associated with the rate constants.
The random walk Metropolis-Hastings took many samples from the posterior distribution of the rate constants. It is, therefore, is straightforward to find the real correlation between each pair of parameters as it is shown in Figure 3. Except for the p4-p5 pair, the correlation between other rate constants is relatively low. The high correlation factor shows that the current kinetic model cannot be fully identified, and the values of p4 and p5 in Table 2 do not correctly represent the rate constants  . From chemistry point, the high correlation coefficient between p4-p5 makes sense because both parameters are belonging to one reaction pathway with equilibrium conditions. In fact, in the equilibrium conditions, the ratio of the rate constants is important. This is why the rate constants of p4 and p5 cannot be identified. This finding is also observed by Yermakova et al.  and they suggest for a reduction in the kinetic model by removing the fifth path.
Following the method presented in  , to show the level of uncertainty associated with the kinetic model of the thermal isomerization of α-pinene, a Monte
Figure 3. Correlation matrix of the thermal isomerization of α-pinene rate constants.
Carlo simulations performed on the posterior samples. 1000 parameters set were randomly selected from the pool of accepted and the kinetic model of Equation (7) was solved for each parameter set. The resulting 95% confidence interval for each component at each experimental time interval (Table 1, first column) is obtained where is shown in Figure 4. The value of ESD has no effect on the model output results as it is shown by blue region in Figure 4 that only includes the uncertainty associated with the parameter values. The blue floating bars show very narrow region for α-pinene (C1) and dipentene (C2) that is consistent with the low standard deviation of 0.039 and 0.037 for p1 and p2, respectively. To take into consideration the model structural error, the ESD of each sample is added to the model output showing the uncertainty associated with the model errors. Therefore, in Figure 4 the red bars are showing the combined uncertainty of parameters as well as model structural.
As it is seen in Figure 4, the duplicate experimental points are used in the parameter estimation without taking their average to take into consideration the real measurement errors. In this way, the experimental data provide better information about the parameter values.
In this study, Bayesian inference and random walk MCMC are applied to estimate the confidence interval of the rate constants of a kinetic model. Thermal isomerization of α-pinene is considered as a case study and the performance of the proposed model is evaluated by obtaining the posterior distribution of its five rate constants plus the error standard deviation of the proposed kinetic model. The results showed that:
1) The Bayesian inference solved by the random walk Metropolis-Hastings technique can successfully estimate the posterior distribution of the rate constants for a proposed kinetic model while experimental data are given.
Figure 4. The thermal isomerization of α-pinene model output confidence interval by considering the parameter and model structural uncertainties.
2) The 95% confidence interval and standard deviation for each rate constants can be obtained from the posterior distribution and its range represents the uncertainty associated with that parameter.
3) The error standard deviation in the likelihood function treated as an unknown parameter in the stochastic model identification approach and its numerical value literally represents the quantified uncertainty of the model.
4) Correlation between each pair of parameters is considered in the MCMC samples. The results of correlation matrix can provide useful information about the identifiability of the proposed kinetic model.
It can be concluded that the confidence interval of rate constants obtained by applying the Bayesian inference concept is an acceptable representation of uncertainties associated with the kinetic model.