Development of a robust nuclear data adjustment method to outliers

. We developed a new nuclear data adjustment method for experimental data containing outliers. This method mitigates the effect of outliers by applying M-estimation, a type of robust estimation, to the conventional nuclear data adjustment method using sensitivity coefficients. Based on the M-estimation, we derived a weighted nuclear data adjustment formula and developed a weight calculation method. The weighted nuclear data adjustment formula was derived by weighting the function to take the extremum of the conventional nuclear data adjustment. The weighting of each nuclear characteristic is calculated from the difference between the measured and calculated values of the nuclear characteristic. This weight calculation method can evaluate the validity of each nuclear characteristic by considering correlations between nuclear characteristics using singular value decomposition. The proposed method and the conventional method were compared and verified by twin experiments. In the twin experiments, the nuclear data were adjusted using experimental data that intentionally included outliers. As a result of twin experiments, it was confirmed that the nuclear data were adjusted robustly and appropriately even with the experimental data containing outliers.


Introduction
Nuclear data adjustment is a technique to reduce the uncertainty of nuclear data by unifying the measured nuclear characteristics in integral experiments and the numerically calculated nuclear characteristics using nuclear data. In conventional nuclear data adjustment methods, the experimental data (measured nuclear characteristics and their uncertainties) are assumed to be appropriate. Note, however, that the adjustment will not be appropriate if experimental data contains outliers because adjusted nuclear data are significantly affected by experimental data. To avoid this issue, the analyst removes experimental data containing outliers based on experience and engineering judgment in the practice of conventional nuclear data adjustment.
For example, in the previous study by Takeda et al. [1], removed data were determined from the difference between measured nuclear characteristic and calculated nuclear characteristic ( / value). Specifically, the authors considered three cases in which nuclear characteristics were removed where the difference exceeded one, two, or three times the total uncertainty (experimental error + analytical uncertainty + nuclear data-induced uncertainty).
As another example, in the adjusted nuclear data ADJ2017 [2], experimental data were removed based on the following points: (i) The data where the difference between measured and calculated values exceeds twice the total uncertainty (2 ); * Corresponding author : y-fukui@fermi.energy.nagoya-u.ac.jp (ii) The data where experimental analysis is considered inappropriate; (iii) The data where the amount of adjustment of nuclear data exceeds its uncertainty (1 ). In particular, the point (iii) poses an inverse problem that searches for the causal experimental data from the adjusted nuclear data, thus it is not easy to specify the experimental data to be removed. In order to accurately specify the group of experimental data that is the cause of the problem, the number of analyses would be very large because inclusion or exclusion of each experimental data will be necessary. Therefore, based on the experience and engineering judgment of the analyst, inappropriate experimental data groups are excluded.
To address the above issues, this study aims to apply a robust estimation method to the nuclear data adjustment. The least median of squares method, the random sample consensus, and the M-estimation are examples of robust estimation methods. In this study, we focus on the M-estimation due to its applicability to Bayesian inference. By introducing the M-estimation to the conventional nuclear data adjustment methods, we develop a new nuclear data adjustment method that is robust to experimental data including outliers.

GLLS
This section describes a conventional nuclear data adjustment method called Generalized Linear Least Square (GLLS), which uses sensitivity coefficients. This method uses the following assumptions:  The perturbation of nuclear characteristics is linear with respect to the perturbation of nuclear data.  The probability distribution of nuclear data and nuclear characteristics follows a multivariate normal distribution. From Bayes' theorem shown in equation (1), the nuclear data adjustment equation is derived [3].
where ⃗ is the nuclear data; ⃗ is the nuclear characteristic measurement value; ( ⃗ ) is the probability distribution (prior probability) of the nuclear data ⃗ ; ( ⃗ ) is the probability distribution of the nuclear characteristic measurement value ⃗ ; ( ⃗ | ⃗ ) is the probability distribution (posterior probability) of the nuclear data when the nuclear characteristic is ⃗ ; ( ⃗ | ⃗ ) is the probability distribution (likelihood) of the nuclear characteristic measurement value when the nuclear data is ⃗ . With the assumption of a multivariate normal distribution for the nuclear data and nuclear characteristics, equations (2)-(4) are derived for ( ⃗ ), ( ⃗ ), and ( ⃗ | ⃗ ), respectively: where ⃗ (0) ∈ ℝ is the prior mean of the nuclear data; (0) ∈ ℝ × is the prior covariance matrix of the nuclear data; is the total number of nuclear data; ⃗ ∈ ℝ is the true value of nuclear characteristics; ∈ ℝ × is the measured covariance matrix; ( ⃗ ) ∈ ℝ is the calculation value of nuclear characteristics for nuclear data ⃗ ; + ∈ ℝ × is the covariance matrix of nuclear characteristic measurement and calculation; is the total number of nuclear characteristics. The posterior probability ( ⃗ | ⃗ ) is obtained by substituting equations (2)-(4) into equation where part of the posterior probability with regard to the nuclear data ⃗ is defined as ( ⃗ ).
Because of the assumption of linearity between nuclear data and nuclear characteristics, the calculated nuclear characteristic value ( ⃗ ) for the nuclear data ⃗ can be calculated using the sensitivity coefficient matrix ∈ ℝ × as follows: The nuclear data ⃗ that maximizes the posterior probability ( ⃗ | ⃗ ) is the "adjusted nuclear data ⃗ (adj) " estimated by Bayes' theorem. To obtain the nuclear data ⃗ that maximizes the posterior probability, we consider the derivative of the function part ( ⃗ ) regard to ⃗ in the posterior probability ( ⃗ | ⃗ ) : The nuclear data ⃗ satisfying equation (9) is the adjusted nuclear data ⃗ (adj) . In addition, by calculating the variance of ⃗ (adj) − ⃗ (0) , the covariance matrix (adj) of adjusted nuclear data can be obtained. From the above, the nuclear data adjustment formulae are obtained as follows:

M-estimation
M-estimation [4] is one of the estimation methods that are robust to outliers. In M-estimation, the regression model is updated by setting the weight ( ) of the observed values from the residual between the observed values and the regression model. The initial regression model (before updating the weight) is determined in advance using the least squares method or other methods. In this method, the influence of outliers on the estimation results can be reduced because the weighting ( ) on the values will be small when the residuals are large. The function that calculates the weight from the residual is called the loss function. An example of the loss function is shown in equation (13).
where is the tolerable error. The estimation using this loss function is called Tukey's Biweight method. [4] In this section, we consider the regression problem for a linear function as an example of M-estimation. Let us consider a problem estimating the parameters and in a regression model (model) = + using pairs of given input and observed output data { , (obs) } (1 ≤ ≤ ).
The ordinary least squares method sets up the simultaneous equations shown in equation (14) and calculates the parameters that minimize the sum of squares of the residuals between the model and the observation: We describe equation (14) in matrix form. If the observed data (obs) , the regression parameter = [ ] T , and the input are used, the following can be obtained for the ordinary least squares method using the pseudo-inverse matrix: where † is the pseudo-inverse of the matrix .
On the other hand, in M-estimation, we set up a simultaneous equation obtained by multiplying the weights . The function to be minimized in this case is shown in equations (20) and (21). The weight of the th data (obs) is calculated by the residual between the calculated value and the observed value by the temporarily set regression model: Here, we describe equation (20) in matrix form using the weight matrix . The weight matrix is a matrix that has the weight in the diagonal: The procedure for applying M-estimation to the problem of regression to a linear function is summarized below.
(i) Initial parameters ( , ) are tentatively determined using the least squares method, etc. (ii) The weights are calculated from the regression model with the initial parameters and the observed values. (iii) The regression parameters are updated using the weights. The procedure for applying M-estimation to nuclear data adjustment is as follows.
(i) The initial parameters are the nuclear data mean (evaluated nuclear data). (ii) The weights are obtained from the calculated and measured nuclear characteristics for the nuclear data mean. (iii) The nuclear data are adjusted using weights. The next section describes the application of Mestimation to nuclear data adjustment in more detail.

Fundamental concept
We consider a nuclear data adjustment method that introduces M-estimation. The basic concept of the proposed method is the adjustment of nuclear data with a focus on nuclear characteristics where the difference between calculated and measured values is within an acceptable uncertainty range. The key points to be developed are described below.
(i) We develop a method to calculate the weighting from the difference between calculated and measured values, taking into account correlations between nuclear characteristics. Typical M-estimation cannot handle correlations between nuclear characteristics because the weights are calculated based on the assumption that the observed data are independent. To solve this problem, we propose a weights calculation method using the singular value decomposition. From the singular value decomposition of the total nuclear characteristic covariance matrix (i.e., characterizing the correlation direction is obtained. This operation allows nuclear characteristics distributions, including correlations, to be treated as independent distributions. The weights can be calculated by evaluating the difference between the calculated and measured values in each base direction. As an example, in the case of the two nuclear characteristics shown in Fig.1, the weighting is larger for the direction of  direction while that is within the scatter for the → 1 direction.

Fig.1 The concept of weights calculation
(ii) We develop a weighted nuclear data adjustment equation for each nuclear characteristics incorporating M-estimation.

Weights calculation
This section describes the setting of the weight matrix , in which the correlations between the nuclear characteristics are taken into account. In this method, the direction (i.e., base) of the correlation is determined from the singular value decomposition of the covariance matrix of the nuclear characteristics, and the weight matrix ( ) is calculated based on the loss function in the obtained base direction appeared in the base matrix . The detailed calculation procedure is shown below.
(2) Obtain the total covariance matrix of the total uncertainty from the sum of the covariance matrices: (3) We obtain the base matrix of the nuclear characteristic correlations by the singular value decomposition. The covariance matrix is decomposed into the orthonormal and the diagonal matrix : where the column elements (5) Calculate the weight matrix ( ) (diagonal matrix) by evaluating the difference between the transformed measured value ⃗ ( ) and the transformed calculated value ( ) using the loss function : where is an arbitrary tuning parameter. In this paper, we set = 1.0. When calculating the inverse of the weight matrix ( )−1 , the diagonal components of ( ) should not contain small values very close to zero. If the diagonal components of ( ) contain very small values, will be a singular matrix and its inverse −1 may not be properly calculated.

Weighted nuclear data adjustment
Like M estimation, the nuclear data adjustment formula is derived by assuming a weighted model equation. As an alternative of the function ( ⃗ ) , which is the minimization target in the conventional nuclear data adjustment, let us consider the function ( ⃗ ), in which nuclear characteristics ⃗ and ( ⃗ ) are weighted: Here, we find that ( ⃗ ) is obtained by replacing the nuclear characteristic uncertainty in the conventional nuclear data adjustment formula to From this result, the conventional nuclear data adjustment formula (10)-(12) is rewritten as: In this paper, we call this method MGLLS (weighted nuclear data adjustment with M estimation incorporated into GLLS).

Verification conditions
We verify the proposed method using a twin experiment, which is one of the verification methods for data assimilation. The procedure for verification in this study is described below.
(2) Numerically calculate the hypothetical measured nuclear characteristics ⃗ using the nuclear data is ⃗ (ref) . (3) Intentionally change a part of the measured nuclear characteristics ⃗ to virtually generate an abnormal value (i.e., outlier). (4) Adjust the nuclear data ⃗ (0) using the measured nuclear characteristics ⃗ , and obtain adjusted nuclear data ⃗ (adj) and its covariance (adj) .
The verification conditions are summarized in Table  1. We analyzed four nuclear characteristics for each of the three experiments. Therefore, the total number of nuclear characteristics is twelve.

Experiments
Pu-Met-Fast-001 (Jezebel) [5] Pu-Met-Fast-006 (Flattop) [5] Pu-Sol-Therm-001 [  To verify the robustness to outliers, we intentionally replaced measured nuclear characteristics ⃗ to an abnormal value. In this study, abnormal experimental data were obtained by changing the sign of the relative difference ∆= ( − )/ for one of the nuclear characteristics, here is a calculated nuclear characteristic for nuclear data ⃗ (0) . In other words, we reversed the positives and negatives of the relative difference. As an example, let us consider the abnormal measured value ′ against the measured value and the calculated value . Because of the positive and negative reversal of its relative perturbations, we can calculate the following equation for the abnormal measured value ′ , In this study, we chose C49/F41 at PMF001 as an outlier because this nuclear characteristic is highly sensitive to the 239 Pu ( , ) cross-section.

Result of adjustment
For experimental data containing outliers, the nuclear data were adjusted using the conventional method (GLLS) and the proposed method (MGLLS). We compare the adjusted nuclear data ⃗ (adj) and its uncertainty after adjustment (i.e., square root of diagonal elements of (adj) ) to reference values ⃗ (ref) .
The cases with different ⃗ (ref) s (#1-#9) are shown in Fig. 2-Fig. 10, respectively. In these figures, the red single dot chain line means the reference nuclear data, i.e., virtual true nuclear data in the twin experiments. When the reference values (red single dot chain line) are within the range of estimation uncertainty (light blue dashed lines), the adjustment can be considered appropriate.
In cases #1,#3-#9, the difference between the reference and adjusted values of nuclear data in GLLS is much larger than its standard deviation. In contrast, in the MGLLS, the difference between the reference and adjusted values of the nuclear data is more consistent with the range of uncertainty (i.e., the standard deviation ).
The result of case #2 shows no significant difference between GLLS and MGLLS. This is because, in case #2, the perturbations of PMF-001-C49/F41 to outliers are small. In other words, since the sign of a small perturbation amount is reversed to generate the outlier of case #2, there is a small difference between the original and the perturbed (including outliers) nuclear characteristics.

Quantitative performance evaluation of the adjustment methods
The adequacy of the adjusted nuclear data is quantitatively evaluated using the 2 value as an index. If the difference between the adjusted value and the reference value is within the uncertainty range of the adjustment, the order of 2 is comparable to that of the degrees of freedom, and the adequacy of the adjustment (or estimation) can be evaluated. The calculation method of the 2 value is as follows: If (adj) is a singular matrix, its inverse matrix cannot be properly calculated. Therefore, a low-rank approximation using singular value decomposition for (adj) is used: where r is the order after the application of the low-rank approximation. And where [: , : ] is slicing to extract the first to r-th column element of matrix , and [: , : ] is the slicing to extract the first to r-th square matrix element of matrix . Using the above singular value decomposition for (adj) , we approximate the 2 values using the following equation: For both adjustment methods, the evaluated 2 values are shown in Fig. 11. The blue circled and the orange x-shaped markers indicate 2 values for the conventional method (GLLS) and the proposed method (MGLLS), respectively. When the adjustment is appropriate, 2 takes the comparable value to the degree of freedom. Because the 239 Pu ( , γ) cross-section is adjusted by 56 group energy bins, the number of degrees of freedom is less or equal to 56 in the present calculation condition. In Fig. 11, the degrees of freedom are indicated by a black dashed line.

Fig. 11 Evaluated value in each case
In cases #1, #3-#9, the 2 values in the GLLS case are extremely large, ranging from 1000 to 100000. On the other hand, the 2 values are comparable to the degrees of freedom (number of nuclear data) in the MGLLS case. The result demonstrates that the estimation results have improved using MGLLS. The reason for almost the same results using GLLS and MGGLS in case #2 is, as mentioned earlier, the small amount of perturbation to generate the outlier.

Conclusions
By incorporating M-estimation into the conventional nuclear data adjustment method, we proposed a robust nuclear data adjustment method for outliers, i.e., MGLLS. In this method, the "adequacy" of each nuclear characteristic is determined from the difference between the measured nuclear characteristic obtained by critical experiments and the calculated nuclear characteristic using nuclear data. Based on the "adequacy," each nuclear characteristic is weighted and the nuclear data is adjusted. Correlations between nuclear characteristics are considered in evaluating this "adequacy." The proposed method was verified through twin experiments on critical experiments including fast and thermal systems. We quantitatively confirmed that the proposed method can robustly and appropriately adjust nuclear data even for experimental data containing inappropriate values.
A future study is an extension of the present method to the Bayesian Monte Carlo (BMC) method [11]. The BMC method uses the nuclear characteristics covariance matrix in the form of an inverse matrix when adjusting nuclear data. Therefore, the weight matrix can be used directly instead of the inverse matrix, thus we can avoid the instability of the inverse matrix calculation of the weight matrix.