# Comprehensive Uncertainty Quantification in Nuclear Safeguards

Authors: E. Bonner; T. Burr; T. Krieger; K. Martin; C. Norman
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

Nuclear safeguards aim to confirm that nuclear materials and activities are used for peaceful purposes. To ensure that States are honoring their safeguards obligations, quantitative conclusions regarding nuclear material inventories and transfers are needed. Statistical analyses used to support these conclusions require uncertainty quantification (UQ), usually by estimating the relative standard deviation (RSD) in random and systematic errors associated with each measurement method. This paper has two main components. First, it reviews why UQ is needed in nuclear safeguards and examines recent efforts to improve both top-down (empirical) UQ and bottom-up (first-principles) UQ for calibration data. Second, simulation is used to evaluate the impact of uncertainty in measurement error RSDs on estimated nuclear material loss detection probabilities in sequences of measured material balances.


---

## Body

## 1. Introduction

Nuclear material accounting (NMA) provides a quantitative basis to detect nuclear material loss or diversion at declared nuclear facilities. NMA involves periodically measuring facility input transfersTin, output transfers Tout, and physical inventory I to compute a material balance (MB) defined for balance period j as MBj=(Ij-1+Tin,j-Tout,j)-Ij, where (Ij-1+Tin,j-Tout,j) is the book inventory. In NMA, one MB or a collection of MBs are tested for the presence of any statistically significant large differences and/or for trends, while allowing for random and systematic errors in variance propagation to estimate the measurement error standard deviation of MBj, σMBj. Similarly, in verification activities done by an inspector, paired operator and inspector data are tested for any large differences and/or for trends [1, 2]. Therefore, both material balance evaluation and verification activities require statistical analyses, which require UQ.In metrology for nuclear safeguards, the term “uncertainty” characterizes the dispersion of estimates of a quantity known as the measurand, which is typically the amount of NM (such as U or Pu) in an item. To measure the amount of NM, both destructive analysis (DA, a sample of item material is analyzed by mass spectrometry in an analytical chemistry laboratory) and nondestructive assay (NDA, an item is assayed by using a neutron or gamma detector) are used. NDA uses calibration and modeling to infer NM mass on the basis of radiation particles, such as neutrons and gammas emitted by the item and registered by detectors. For any measurement technique, one can use a first-principles physics-based or “bottom-up” approach to UQ by considering each key step and assumption of the particular method. Alternatively, one can take an empirical or “top-down” approach to UQ, for example, by comparing assay results on the same or similar items by multiple laboratories and/or calibration periods.A well-known guide for bottom-up UQ is the Guide to the Expression of Uncertainty in Measurement (GUM, [3]). The GUM also briefly mentions top-down UQ in the context of applying analysis of variance (ANOVA, [4]) to data from interlaboratory studies. Although the GUM is useful, it is being revised because it has known limitations [5–7]. For example, the GUM provides little technical guidance regarding calibration as a type of bottom-up UQ or regarding top-down UQ [5–8]. The GUM also mixes Bayesian with non-Bayesian concepts. In a Bayesian approach, all quantities, including the true measurand value, are regarded as random. In a non-Bayesian (frequentist) approach, some quantities are regarded as random and other quantities, such as the true value of the measurand, are regarded as unknown constants. This paper uses both Bayesian and non-Bayesian concepts but specifies when each is in effect. For example, in the Bayesian approach to top-down UQ in Section 3, the true RSD values are regarded as being random.In NDA safeguards applications, the facility operator declares the NM mass of each item. Then, some of those items are randomly selected for NDA verification measurement by inspectors. This is a challenging NDA application because often the detector is brought to the facility where ambient conditions can vary over time and because the items to be assayed are often heterogeneous in some way and/or are different from the items that were used to calibrate/validate and assess uncertainty in the NDA method. Because of such challenges, “dark uncertainty” [9] can be large, as is evident whenever bottom-up UQ predicts smaller measurement error RSDs than are observed in top-down UQ [1]. The RSD of an assay method is often defined as the reproducibility standard deviation as estimated in an interlaboratory comparison. As shown in Section 3, comparing NDA verification measurements to the operator’s DA measurements can be regarded as a special case of an interlaboratory evaluation [10–12].For top-down UQ applied to NM measurements of the same item by both the operator (often using DA) and the inspector (often using NDA), this paper describes an existing and a new approach to separately estimate operator and inspector systematic and random error variance components. Systematic and random error components must be separated because their modes of propagation are different (Section4). Currently, random error variance estimates (from paired data) are based on Grubbs’ estimator or variations of Grubbs’ estimator, which was originally developed by Grubbs to estimate random error variance separately for each of the two methods applied to each of several items, without repetition of measurement by either method [13, 14]. In Section 3, Grubbs’ estimator, constrained versions of Grubbs’ estimator, and a Bayesian alternative [7] are described; the Bayesian option easily allows for parameter constraints and prior information regarding the relative magnitudes of variance components to be exploited to improve top-down UQ.This paper is organized as follows. Section2 provides a background on bottom-up UQ for NDA, describes a gamma-based NDA example and a neutron-based NDA example, and illustrates why simulation is necessary for improved UQ for calibration data. Section 3 reviews currently used top-down UQ and describes a new Bayesian option [7] that applies approximate Bayesian computation. Section 4 provides a new simulation study assessing the sensitivity of estimated NM loss detection probabilities to estimation errors in measurement error RSDs. Section 5 concludes with a summary.


## 2. Bottom-Up UQ

For bottom-up UQ, the GUM [3] assumes that the measured value can be expressed using a measurand equation that relates input quantities (data collected during the measurement process and relevant fundamental nuclear data such as attenuation corrections) to the output (the final measurement value). The GUM’s main technical tool is a first-order Taylor approximation applied to the measurand equation (1)Y=fX1,X2,…,XN,which relates input quantities X1,X2,…,XN (regarded as random) to the measurandY (also regarded as random). The input quantities can include estimates of other measurands or of calibration parameters, so (1) is quite general. The variance of each X and σi2 and any covariances, σiσjρi,j, between pairs of X’s are then propagated using the Taylor approximation to obtain σY2≈∑i=1N∂f/∂xi2σi2+2∑i=1N-1∑j=i+1N(∂f/∂xi)(∂f/∂xj)σiσjρi,j (or using simulation if the Taylor approximation is not sufficiently accurate) to estimate the variance in Y, σY2.According to (1), the estimated value Y of the measurand is a random variable, regardless of whether the left side of (1) is expressed asY (as in a typical Bayesian approach) or as μ^Y (as in a typical non-Bayesian setting). The hat notation is a frequentist convention for denoting an estimator, so μ^Y is an estimate of μY, and μY (which is also denoted as yT, where “T” denotes the true value) denotes the unknown true value of the measurand.


### 2.1. Calibration UQ as an Example of Bottom-Up UQ

Typically, calibration is performed using reference materials having nominal measurand values (known to within a relatively small uncertainty), and then, in the case of linear calibration, (1) can be reexpressed as μ^Y=β^O+β^1X, where μ^Y is the estimated measurand value, β^0 and β^1 are parameters estimated from calibration data, X is the net count rate (usually the net gamma or net neutron count rate in NDA; see Section 2.3), and the three inputs in mapping to (1) are X1=β^0, X2=β^1, and X3=X. The estimates β^0 and β^1 will vary in predictable ways (Sections 2.2 and 2.3) across repeats of the calibration.The convention in statistical literature to reverse the roles ofX and Y from that in GUM’s equation (1) will be followed here, so X denotes the quantity to be inferred (the measurand value) and Y denotes the detected net radiation count rate. Then, in the case of reverse regression (see below), (1) can be expressed as (2)X=gY1,Y2,…,YN=α^0+α^1Y,identifying Y1=β^0, Y2=β^1, and Y3=Y. Following calibration on data consisting of n(xi,yi) pairs (lowercase denotes observed value of a random variable), the three “input quantities” Y1=β^0, Y2=β^1, and Y3=Y have variances and covariances that can be estimated. However, in most applications of calibration in NDA, accurate estimation of these variances and covariances requires simulation because analytical approximations as described in Section 2.2 have been shown to be inadequate (see Section 2.3).Expressing (2) as X=g(Y1,Y2,…,YN)=α^0+α^1Y indicates how the estimate X is computed and how to assign systematic and random error variances to X. For example, and to introduce notation used in top-down UQ (Section 3), one could express the estimate as X=μX+S+R, where μX denotes the true value of the measurand,S denotes systematic error due to estimation error in the fitted slope and intercept and/or due to correlations among the inputs, and R denotes random error. If there are no correlations among the inputs but only estimation error in the fitted slope and intercept during calibration, then expressions for the variances of S andR, denoted as σS2 and σR2, respectively, can be given (Sections 2.3 and 3), which allow comparison between bottom-up UQ and top-down UQ. The GUM does not discuss calibration in much detail; instead, the GUM applies propagation of variance to the steps modeled in (1), which sometimes leads to a defensible estimate of the combined variances of S and R. The GUM does not attempt to separately estimate the variances of S andR, but such separation is needed in some applications, such as assigning an uncertainty to a sum of measurand estimates ([15] and Section 4).


### 2.2. Extension of Standard Regression Results to Calibration

One way to express that the net count rate depends on the true measurand value is(3)Y=β0+β1XTrue+RY,which is a typical model used in regression when there is negligible error in X. If errors in predictors cannot be ignored, (3) should be modified; however, one can still regress measured Y on measuredX, so, in effect, (3) can be reexpressed as Y=β~0+β~1X+R~Y, where the tildes denote that parameter values and the random error are different from those in (3).In inverse calibration, (3) is used, and one inverts the fitted model using β^0 and β^1 to use future measured ytest to predict enrichment using(4)x^test=ytest-β^0β^1,which is regression followed by inversion. An alternative model to (3) is reverse calibration:(5)X=α0+α1YTrue+RX.In reverse calibration, (5) expresses the measurand X as a function of the true net count rate yT, but in practice one must regress X on Y=Ytrue+RY, where RY is a random error. As an aside, this paper does not consider models with systematic errors such as X=XTrue+RX+SX or Y=Ytrue+RY+SY. Cheng and Van Ness [16] point out that any additive systematic errors in X or Y could be absorbed into β0 and α0, respectively; however, any systematic errors in the X values used for calibration would remain a part of the total uncertainty.Both inverse and reverse calibrations involve ratios of random variables, which can be problematic [7, 17]. In inverse calibration, the solution in (4) involves division by the random variable β^1, which has a normal distribution under typical modeling assumptions. Williams [18] notes that x^test in (4) has infinite variance even if the expected value of β^1 is nonzero, due to division by a normal random variable [19], and hence has infinite mean squared error, while the reverse estimator has finite variance and mean squared error. In reverse calibration, the least squares solution x^T=α^y=YcalTXcal/YcalTYcaly also involves division of random variables (Xcal is the vector or matrix of X values used in calibration and Ycal is the vector of Y values in calibration). Experience suggests that one can develop adequate approximations for the ratio of random variables when the ratio is almost certain to be far from infinity or zero [19]. Ignoring errors in predictors, [20] uses the following common approximation for the variance of the ratio of random variables U andV:(6)varUV≈EU2EV2varUEU2+varVEV2-2covU,VEUEV,where E denotes expectation to derive the approximation (for inverse calibration) for variance due to uncertainty in the estimated calibration coefficients β^0 and β^1 and in the test measurement Ytest: (7)varX^≈1β12σRY2n+xtest-x-2σRY2∑i=1nx~i2,where x~i=xi-x- and x- is the mean of the x values in the calibration data. To apply (7), β12 and σRY2 are estimated from the calibration data (assuming (3) in forward calibration or the alternate version of (3), Y=β~0+β~1X+R~Y). Equation (7) is almost the same as the corresponding well-known result for regression; the only differences are the swapping of the roles for x andy and the appearance of β12 in the denominator. For reverse regression, [20] derives varx^≈((∑i=1nx~i2)/n-2)1/n+ytest-y-2/∑i=1ny~i2, where y~i=yi-y- and x~i=xi-x-. Reference [20] also showed the long-term bias Binverse≈(x-x-)σRY2/β12Sxx for inverse calibration and Breverse≈-(x-x-)/(1+β12Sxx/(n-1)σRY2) for reverse calibration, where Sxx=∑i=1n(xi-x-)2. Notice that Binverse decreases as n increases (because Sxx increase as n increases), but Breverse does not decrease as n increases; however, recall that, in NDA applications,n is small, usually 3 to 10.A common summary performance measure of an estimator combines squared bias and repeatability variance defined as RMSE = (repeatability variance) + (bias)2; that is, RMSE=E{X^-Xtrue}2=E{X^-EX^}2+{EX^-Xtrue}2, where E denotes the expected value (i.e., the first moment of the underlying probability distribution) [7]. Some technical details arise regarding the best model fitting approach if the predictor Y is measured with nonnegligible error. In addition, there is controversy regarding the relative merits of inverse and reverse calibration [7, 17, 21, 22]. Simulation can be used to choose between inverse and reverse calibration, because simulation provides accurate UQ (such as RMSE estimation) for both options. In simulations for NDA calibration, errors in the standard reference materials’ nominal values (X’s) are usually small compared to errors in the instrument responsesY’s, which are possibly adjusted by using adjustment factors that have uncertainty (see Section 2.3).


### 2.3. Summary of Recent NDA Examples

Recent publications have used simulation to assess the adequacy of (7) in the context of NM measurements by gamma detection [7, 17] and neutron detection [1, 7, 23, 24].


#### 2.3.1. Enrichment Meter Principle (EMP)

The EMP aims to infer the fraction of235U in U (enrichment, defined as atom percent of 235U in an item) by measuring the count rate of the strongest-intensity direct (full-energy) gamma from decay of 235U, which is emitted at 185.7 keV [7, 25, 26]. The EMP makes three key assumptions: (1) the detector field of view into each item is the same as that in the calibration items, (2) the item is homogeneous with respect to both 235U enrichment and chemical composition, and (3) the container attenuation of gamma-rays is the same as or similar to that in the calibration items, so empirical correction factors have modest impact and are reasonably effective. If these three assumptions are approximately met, the enrichment of 235U in the U is directly proportional to the count rate of the 185.7 keV gamma-rays emitted from the item. It has been shown empirically that, under good measurement conditions, the EMP can have a random error RSD of less than 0.5% and long-term bias of less than 1% relative to the true value, depending on the specific implementation of the EMP. Implementation details include features such as the detector resolution, stability, and extent of corrections needed to adjust items to calibration conditions. However, in some EMP applications, the random error RSD can be larger than bottom-up UQ predicts and larger than the 0.5% goal. For example, assay of the 235U mass in a stratum of UO2 drums suggests that there is a larger-than-anticipated random RSD [17].


#### 2.3.2. Uranium Neutron Coincidence Collar (UNCL)

The UNCL uses an active neutron source to induce fission in235U in fresh fuel assemblies [27]. Neutrons from fission are emitted in short bursts of time and so exhibit non-Poisson bursts in detected count rates. Neutron coincidence counting is used to measure the “doubles” neutron coincident rateY, which can be used to estimate the linear density of 235U in a fuel assembly (g-235U/cm) using calibration parameters, a1 and a2. The coincident rate Y is the observed rate of observing two neutrons in very short time gates, each of approximately 10-6sec, and is attributable to fission events. The equation commonly used to convert the measured doubles rate Y to an estimate ofX (grams 235U per cm) is X=kY/(a1-a2kY), where a1 and a2 are calibration parameters, and k=k0k1k2k3k4k5 is a product of correction factors that adjust Y to item-, detector-, and source-specific conditions in the calibration [27]. Therefore, X=kY/(a1-a2kY) is a special case of GUM’s equation (1) (with X and Y reversed), where the two calibration parameters a1 and a2 and the 6 correction factors k0,k1,k2,k3,k4, and k5 are among X’s in (1).Reference [23] showed that calibration is most effective (leading to smallest RMSE in X^) if there is no adjustment for errors in the predictor kY and that errors in k0,k1,k2,k3,k4, and k5, in k=k0k1k2k3k4k5, should be included in synthetic calibration data. Note that, by working with 1/X and 1/Y, one can convert X=kY/(a1-a2kY) to one that is linear in the transformed predictor 1/Y.


#### 2.3.3. Main Results for Sections2.3.1 and 2.3.2

The main results for Sections2.3.1 and 2.3.2 can be summarized in four main points as follows.(1) If possible, both classical (see (2)) and reverse (see (3)) regression methods should be compared; however, reverse regression tends to do either as well as or better than classical regression. Analytical approximations such as (7) have been shown not to be sufficiently adequate in some settings, so simulation is recommended to compare classical and reverse regression and to estimate variance components in X=μX+S+R (Section 3).(2) Error sources that are expected to be present in test measurements, such as container thickness measurements, can be simulated in synthetic calibration data. Such error sources often lead to item-specific biases (Burr et al., 2016).(3) If reverse regression is used, then there is no need to adjust for errors in the predictors Y in (3). If inverse regression is used, then it is better to adjust for errors in predictors.(4) Figure 1 plots (a) the observed and predicted bias and (b) the observed and predicted RMSE in a generic NDA example involving either gamma or neutron counting. It is not well known that calibration applications lead to bias, and [7, 17] showed that the bias cannot be easily removed, because measurement errors obscure the true measurand value and hence the true bias. Note in Figure 1(a) that the observed bias (in simulated data) is not in close agreement with the predicted bias, which is obtained from the expressions in Section 2.2. Therefore, long-term bias should be estimated using simulation rather than relying on the approximate expressions Binverse≈(x-x-)σRY2/β12Sxx for inverse calibration and Breverse≈-(x-x-)/(1+β12Sxx/(n-1)σRY2) for reverse calibration, where Sxx=∑i=1n(xi-x-)2. Similarly, Figure 1(b) illustrates that the observed RMSE is not well predicted by the expressions in Section 2.2, so, again, simulation is needed for adequate estimation of the RMSE. Note that the smallest RMSE is for reverse regression. Burr et al. [7, 17] show that reverse regression tends to have smaller RMSE than inverse regression but that if inverse regression is used, then methods to adjust for errors in predictors should be used.Figure 1
Simulation results for both inverse and reverse regression for (a) the observed and predicted bias versusXTrue and (b) the observed and predicted RMSE versus XTrue.
(a)
(b)Figure1 summarizes the results of 105 simulations of Y=β0+β1XTrue+RY=1+0.1XTrue+RY with δX=0.01 and δY=0.15, using 5 (x,y) calibration pairs (with x scaled to lie in (0,1), at 0, 0.25, 0.5, 0.75, and 1) and 10 testing pairs as shown. All simulations in this paper are done in R [25].


## 2.1. Calibration UQ as an Example of Bottom-Up UQ

Typically, calibration is performed using reference materials having nominal measurand values (known to within a relatively small uncertainty), and then, in the case of linear calibration, (1) can be reexpressed as μ^Y=β^O+β^1X, where μ^Y is the estimated measurand value, β^0 and β^1 are parameters estimated from calibration data, X is the net count rate (usually the net gamma or net neutron count rate in NDA; see Section 2.3), and the three inputs in mapping to (1) are X1=β^0, X2=β^1, and X3=X. The estimates β^0 and β^1 will vary in predictable ways (Sections 2.2 and 2.3) across repeats of the calibration.The convention in statistical literature to reverse the roles ofX and Y from that in GUM’s equation (1) will be followed here, so X denotes the quantity to be inferred (the measurand value) and Y denotes the detected net radiation count rate. Then, in the case of reverse regression (see below), (1) can be expressed as (2)X=gY1,Y2,…,YN=α^0+α^1Y,identifying Y1=β^0, Y2=β^1, and Y3=Y. Following calibration on data consisting of n(xi,yi) pairs (lowercase denotes observed value of a random variable), the three “input quantities” Y1=β^0, Y2=β^1, and Y3=Y have variances and covariances that can be estimated. However, in most applications of calibration in NDA, accurate estimation of these variances and covariances requires simulation because analytical approximations as described in Section 2.2 have been shown to be inadequate (see Section 2.3).Expressing (2) as X=g(Y1,Y2,…,YN)=α^0+α^1Y indicates how the estimate X is computed and how to assign systematic and random error variances to X. For example, and to introduce notation used in top-down UQ (Section 3), one could express the estimate as X=μX+S+R, where μX denotes the true value of the measurand,S denotes systematic error due to estimation error in the fitted slope and intercept and/or due to correlations among the inputs, and R denotes random error. If there are no correlations among the inputs but only estimation error in the fitted slope and intercept during calibration, then expressions for the variances of S andR, denoted as σS2 and σR2, respectively, can be given (Sections 2.3 and 3), which allow comparison between bottom-up UQ and top-down UQ. The GUM does not discuss calibration in much detail; instead, the GUM applies propagation of variance to the steps modeled in (1), which sometimes leads to a defensible estimate of the combined variances of S and R. The GUM does not attempt to separately estimate the variances of S andR, but such separation is needed in some applications, such as assigning an uncertainty to a sum of measurand estimates ([15] and Section 4).


## 2.2. Extension of Standard Regression Results to Calibration

One way to express that the net count rate depends on the true measurand value is(3)Y=β0+β1XTrue+RY,which is a typical model used in regression when there is negligible error in X. If errors in predictors cannot be ignored, (3) should be modified; however, one can still regress measured Y on measuredX, so, in effect, (3) can be reexpressed as Y=β~0+β~1X+R~Y, where the tildes denote that parameter values and the random error are different from those in (3).In inverse calibration, (3) is used, and one inverts the fitted model using β^0 and β^1 to use future measured ytest to predict enrichment using(4)x^test=ytest-β^0β^1,which is regression followed by inversion. An alternative model to (3) is reverse calibration:(5)X=α0+α1YTrue+RX.In reverse calibration, (5) expresses the measurand X as a function of the true net count rate yT, but in practice one must regress X on Y=Ytrue+RY, where RY is a random error. As an aside, this paper does not consider models with systematic errors such as X=XTrue+RX+SX or Y=Ytrue+RY+SY. Cheng and Van Ness [16] point out that any additive systematic errors in X or Y could be absorbed into β0 and α0, respectively; however, any systematic errors in the X values used for calibration would remain a part of the total uncertainty.Both inverse and reverse calibrations involve ratios of random variables, which can be problematic [7, 17]. In inverse calibration, the solution in (4) involves division by the random variable β^1, which has a normal distribution under typical modeling assumptions. Williams [18] notes that x^test in (4) has infinite variance even if the expected value of β^1 is nonzero, due to division by a normal random variable [19], and hence has infinite mean squared error, while the reverse estimator has finite variance and mean squared error. In reverse calibration, the least squares solution x^T=α^y=YcalTXcal/YcalTYcaly also involves division of random variables (Xcal is the vector or matrix of X values used in calibration and Ycal is the vector of Y values in calibration). Experience suggests that one can develop adequate approximations for the ratio of random variables when the ratio is almost certain to be far from infinity or zero [19]. Ignoring errors in predictors, [20] uses the following common approximation for the variance of the ratio of random variables U andV:(6)varUV≈EU2EV2varUEU2+varVEV2-2covU,VEUEV,where E denotes expectation to derive the approximation (for inverse calibration) for variance due to uncertainty in the estimated calibration coefficients β^0 and β^1 and in the test measurement Ytest: (7)varX^≈1β12σRY2n+xtest-x-2σRY2∑i=1nx~i2,where x~i=xi-x- and x- is the mean of the x values in the calibration data. To apply (7), β12 and σRY2 are estimated from the calibration data (assuming (3) in forward calibration or the alternate version of (3), Y=β~0+β~1X+R~Y). Equation (7) is almost the same as the corresponding well-known result for regression; the only differences are the swapping of the roles for x andy and the appearance of β12 in the denominator. For reverse regression, [20] derives varx^≈((∑i=1nx~i2)/n-2)1/n+ytest-y-2/∑i=1ny~i2, where y~i=yi-y- and x~i=xi-x-. Reference [20] also showed the long-term bias Binverse≈(x-x-)σRY2/β12Sxx for inverse calibration and Breverse≈-(x-x-)/(1+β12Sxx/(n-1)σRY2) for reverse calibration, where Sxx=∑i=1n(xi-x-)2. Notice that Binverse decreases as n increases (because Sxx increase as n increases), but Breverse does not decrease as n increases; however, recall that, in NDA applications,n is small, usually 3 to 10.A common summary performance measure of an estimator combines squared bias and repeatability variance defined as RMSE = (repeatability variance) + (bias)2; that is, RMSE=E{X^-Xtrue}2=E{X^-EX^}2+{EX^-Xtrue}2, where E denotes the expected value (i.e., the first moment of the underlying probability distribution) [7]. Some technical details arise regarding the best model fitting approach if the predictor Y is measured with nonnegligible error. In addition, there is controversy regarding the relative merits of inverse and reverse calibration [7, 17, 21, 22]. Simulation can be used to choose between inverse and reverse calibration, because simulation provides accurate UQ (such as RMSE estimation) for both options. In simulations for NDA calibration, errors in the standard reference materials’ nominal values (X’s) are usually small compared to errors in the instrument responsesY’s, which are possibly adjusted by using adjustment factors that have uncertainty (see Section 2.3).


## 2.3. Summary of Recent NDA Examples

Recent publications have used simulation to assess the adequacy of (7) in the context of NM measurements by gamma detection [7, 17] and neutron detection [1, 7, 23, 24].


### 2.3.1. Enrichment Meter Principle (EMP)

The EMP aims to infer the fraction of235U in U (enrichment, defined as atom percent of 235U in an item) by measuring the count rate of the strongest-intensity direct (full-energy) gamma from decay of 235U, which is emitted at 185.7 keV [7, 25, 26]. The EMP makes three key assumptions: (1) the detector field of view into each item is the same as that in the calibration items, (2) the item is homogeneous with respect to both 235U enrichment and chemical composition, and (3) the container attenuation of gamma-rays is the same as or similar to that in the calibration items, so empirical correction factors have modest impact and are reasonably effective. If these three assumptions are approximately met, the enrichment of 235U in the U is directly proportional to the count rate of the 185.7 keV gamma-rays emitted from the item. It has been shown empirically that, under good measurement conditions, the EMP can have a random error RSD of less than 0.5% and long-term bias of less than 1% relative to the true value, depending on the specific implementation of the EMP. Implementation details include features such as the detector resolution, stability, and extent of corrections needed to adjust items to calibration conditions. However, in some EMP applications, the random error RSD can be larger than bottom-up UQ predicts and larger than the 0.5% goal. For example, assay of the 235U mass in a stratum of UO2 drums suggests that there is a larger-than-anticipated random RSD [17].


### 2.3.2. Uranium Neutron Coincidence Collar (UNCL)

The UNCL uses an active neutron source to induce fission in235U in fresh fuel assemblies [27]. Neutrons from fission are emitted in short bursts of time and so exhibit non-Poisson bursts in detected count rates. Neutron coincidence counting is used to measure the “doubles” neutron coincident rateY, which can be used to estimate the linear density of 235U in a fuel assembly (g-235U/cm) using calibration parameters, a1 and a2. The coincident rate Y is the observed rate of observing two neutrons in very short time gates, each of approximately 10-6sec, and is attributable to fission events. The equation commonly used to convert the measured doubles rate Y to an estimate ofX (grams 235U per cm) is X=kY/(a1-a2kY), where a1 and a2 are calibration parameters, and k=k0k1k2k3k4k5 is a product of correction factors that adjust Y to item-, detector-, and source-specific conditions in the calibration [27]. Therefore, X=kY/(a1-a2kY) is a special case of GUM’s equation (1) (with X and Y reversed), where the two calibration parameters a1 and a2 and the 6 correction factors k0,k1,k2,k3,k4, and k5 are among X’s in (1).Reference [23] showed that calibration is most effective (leading to smallest RMSE in X^) if there is no adjustment for errors in the predictor kY and that errors in k0,k1,k2,k3,k4, and k5, in k=k0k1k2k3k4k5, should be included in synthetic calibration data. Note that, by working with 1/X and 1/Y, one can convert X=kY/(a1-a2kY) to one that is linear in the transformed predictor 1/Y.


### 2.3.3. Main Results for Sections2.3.1 and 2.3.2

The main results for Sections2.3.1 and 2.3.2 can be summarized in four main points as follows.(1) If possible, both classical (see (2)) and reverse (see (3)) regression methods should be compared; however, reverse regression tends to do either as well as or better than classical regression. Analytical approximations such as (7) have been shown not to be sufficiently adequate in some settings, so simulation is recommended to compare classical and reverse regression and to estimate variance components in X=μX+S+R (Section 3).(2) Error sources that are expected to be present in test measurements, such as container thickness measurements, can be simulated in synthetic calibration data. Such error sources often lead to item-specific biases (Burr et al., 2016).(3) If reverse regression is used, then there is no need to adjust for errors in the predictors Y in (3). If inverse regression is used, then it is better to adjust for errors in predictors.(4) Figure 1 plots (a) the observed and predicted bias and (b) the observed and predicted RMSE in a generic NDA example involving either gamma or neutron counting. It is not well known that calibration applications lead to bias, and [7, 17] showed that the bias cannot be easily removed, because measurement errors obscure the true measurand value and hence the true bias. Note in Figure 1(a) that the observed bias (in simulated data) is not in close agreement with the predicted bias, which is obtained from the expressions in Section 2.2. Therefore, long-term bias should be estimated using simulation rather than relying on the approximate expressions Binverse≈(x-x-)σRY2/β12Sxx for inverse calibration and Breverse≈-(x-x-)/(1+β12Sxx/(n-1)σRY2) for reverse calibration, where Sxx=∑i=1n(xi-x-)2. Similarly, Figure 1(b) illustrates that the observed RMSE is not well predicted by the expressions in Section 2.2, so, again, simulation is needed for adequate estimation of the RMSE. Note that the smallest RMSE is for reverse regression. Burr et al. [7, 17] show that reverse regression tends to have smaller RMSE than inverse regression but that if inverse regression is used, then methods to adjust for errors in predictors should be used.Figure 1
Simulation results for both inverse and reverse regression for (a) the observed and predicted bias versusXTrue and (b) the observed and predicted RMSE versus XTrue.
(a)
(b)Figure1 summarizes the results of 105 simulations of Y=β0+β1XTrue+RY=1+0.1XTrue+RY with δX=0.01 and δY=0.15, using 5 (x,y) calibration pairs (with x scaled to lie in (0,1), at 0, 0.25, 0.5, 0.75, and 1) and 10 testing pairs as shown. All simulations in this paper are done in R [25].


## 2.3.1. Enrichment Meter Principle (EMP)

The EMP aims to infer the fraction of235U in U (enrichment, defined as atom percent of 235U in an item) by measuring the count rate of the strongest-intensity direct (full-energy) gamma from decay of 235U, which is emitted at 185.7 keV [7, 25, 26]. The EMP makes three key assumptions: (1) the detector field of view into each item is the same as that in the calibration items, (2) the item is homogeneous with respect to both 235U enrichment and chemical composition, and (3) the container attenuation of gamma-rays is the same as or similar to that in the calibration items, so empirical correction factors have modest impact and are reasonably effective. If these three assumptions are approximately met, the enrichment of 235U in the U is directly proportional to the count rate of the 185.7 keV gamma-rays emitted from the item. It has been shown empirically that, under good measurement conditions, the EMP can have a random error RSD of less than 0.5% and long-term bias of less than 1% relative to the true value, depending on the specific implementation of the EMP. Implementation details include features such as the detector resolution, stability, and extent of corrections needed to adjust items to calibration conditions. However, in some EMP applications, the random error RSD can be larger than bottom-up UQ predicts and larger than the 0.5% goal. For example, assay of the 235U mass in a stratum of UO2 drums suggests that there is a larger-than-anticipated random RSD [17].


## 2.3.2. Uranium Neutron Coincidence Collar (UNCL)

The UNCL uses an active neutron source to induce fission in235U in fresh fuel assemblies [27]. Neutrons from fission are emitted in short bursts of time and so exhibit non-Poisson bursts in detected count rates. Neutron coincidence counting is used to measure the “doubles” neutron coincident rateY, which can be used to estimate the linear density of 235U in a fuel assembly (g-235U/cm) using calibration parameters, a1 and a2. The coincident rate Y is the observed rate of observing two neutrons in very short time gates, each of approximately 10-6sec, and is attributable to fission events. The equation commonly used to convert the measured doubles rate Y to an estimate ofX (grams 235U per cm) is X=kY/(a1-a2kY), where a1 and a2 are calibration parameters, and k=k0k1k2k3k4k5 is a product of correction factors that adjust Y to item-, detector-, and source-specific conditions in the calibration [27]. Therefore, X=kY/(a1-a2kY) is a special case of GUM’s equation (1) (with X and Y reversed), where the two calibration parameters a1 and a2 and the 6 correction factors k0,k1,k2,k3,k4, and k5 are among X’s in (1).Reference [23] showed that calibration is most effective (leading to smallest RMSE in X^) if there is no adjustment for errors in the predictor kY and that errors in k0,k1,k2,k3,k4, and k5, in k=k0k1k2k3k4k5, should be included in synthetic calibration data. Note that, by working with 1/X and 1/Y, one can convert X=kY/(a1-a2kY) to one that is linear in the transformed predictor 1/Y.


## 2.3.3. Main Results for Sections2.3.1 and 2.3.2

The main results for Sections2.3.1 and 2.3.2 can be summarized in four main points as follows.(1) If possible, both classical (see (2)) and reverse (see (3)) regression methods should be compared; however, reverse regression tends to do either as well as or better than classical regression. Analytical approximations such as (7) have been shown not to be sufficiently adequate in some settings, so simulation is recommended to compare classical and reverse regression and to estimate variance components in X=μX+S+R (Section 3).(2) Error sources that are expected to be present in test measurements, such as container thickness measurements, can be simulated in synthetic calibration data. Such error sources often lead to item-specific biases (Burr et al., 2016).(3) If reverse regression is used, then there is no need to adjust for errors in the predictors Y in (3). If inverse regression is used, then it is better to adjust for errors in predictors.(4) Figure 1 plots (a) the observed and predicted bias and (b) the observed and predicted RMSE in a generic NDA example involving either gamma or neutron counting. It is not well known that calibration applications lead to bias, and [7, 17] showed that the bias cannot be easily removed, because measurement errors obscure the true measurand value and hence the true bias. Note in Figure 1(a) that the observed bias (in simulated data) is not in close agreement with the predicted bias, which is obtained from the expressions in Section 2.2. Therefore, long-term bias should be estimated using simulation rather than relying on the approximate expressions Binverse≈(x-x-)σRY2/β12Sxx for inverse calibration and Breverse≈-(x-x-)/(1+β12Sxx/(n-1)σRY2) for reverse calibration, where Sxx=∑i=1n(xi-x-)2. Similarly, Figure 1(b) illustrates that the observed RMSE is not well predicted by the expressions in Section 2.2, so, again, simulation is needed for adequate estimation of the RMSE. Note that the smallest RMSE is for reverse regression. Burr et al. [7, 17] show that reverse regression tends to have smaller RMSE than inverse regression but that if inverse regression is used, then methods to adjust for errors in predictors should be used.Figure 1
Simulation results for both inverse and reverse regression for (a) the observed and predicted bias versusXTrue and (b) the observed and predicted RMSE versus XTrue.
(a)
(b)Figure1 summarizes the results of 105 simulations of Y=β0+β1XTrue+RY=1+0.1XTrue+RY with δX=0.01 and δY=0.15, using 5 (x,y) calibration pairs (with x scaled to lie in (0,1), at 0, 0.25, 0.5, 0.75, and 1) and 10 testing pairs as shown. All simulations in this paper are done in R [25].


## 3. Top-Down UQ to Estimate Variance Components

In facilities under international safeguards agreements, inspectors measure randomly selected items to monitor for possible data falsification by the operator that could mask NM diversion [1, 28]. These paired (O, I) data (O denotes operator measurement; I denotes inspector measurement) are assessed using one-item-at-a-time testing to detect significant differences and also by using an average of the operator-inspector values to detect trends in the context of material balance evaluation [1]. Conclusions from such an assessment depend on the assumed measurement error model and associated random and systematic uncertainty components, so it is important to perform effective UQ [1, 7, 17, 28].The paired (O, I) data are collected during relatively short (one week) inspections that occur once or a few times per year, and then several years of paired (O,I) inspection data are included in top-down UQ. The measurement error model must account for variation within and between groups, where, in this context, a group is an inspection period. A typical top-down model used for additive errors for the inspector (I) (and similarly for the operator O) is(8)Iij=μij+SIi+RIij,where Iij is the inspector’s measured value of item j (from 1 to n) in group i (from 1 to g), μij is the true but unknown value of item j from group i,RIij~N(0,σRI2) is a random error on item j from group i, and SIi~N(0,σSI2) is a short-term systematic error in group i [28]. The error variance components σSI2 and σRI2 can be estimated using a specialized version of random-effects one-way ANOVA described in Section 3.1. NDA measurements often have larger uncertainty at larger true values, which implies a multiplicative rather than an additive error model. However, provided that the individual RSDs are fairly small, resulting in a total RSD of approximately 10% or less, a multiplicative error model such as Iij=μij(1+SIi+RIij) can be analyzed in the same manner as an additive error model, by analyzing on the log scale [1, 2]. Therefore, for brevity, only an additive error model such as in (8) is presented here. Bonner et al. [1] provide new expressions for a multiplicative error model that should be used if the total RSD is approximately 10% or larger. Note that one could write (8) in more cluttered notation as YIij=μij+SIi+RIij. That is, Iij is the inspector’s measured value of item j, which is obtained using various inputs, denoted with X’s on the right side of (1). And one could also consider other error models, such as error models that allow for nonconstant absolute or relative random and/or systematic SD [28, 29].The GUM [3] briefly describes ANOVA in the context of top-down UQ using measurement results from multiple laboratories and/or assay methods to measure the same measurand; however, the GUM is mostly concerned with bottom-up UQ. The GUM does not explicitly present any measurement error models such as (8) but only considers the model for the measurand, (1). However, the GUM implicitly endorses the notion of a measurement error model (or “observation equation,” [14]) such as (8) in its top-down UQ. Note that if total measurement error is partitioned into random and systematic components, then the variance of a sum of n measured NM values (which is often needed in safeguards assessments; see Section 4) is n2σSI2+nσRI2 for an additive model such as (8).To illustrate, Figure2 plots d=O-I data simulated from (8) with n=10,g = 5,σRO=1,σSO=0.1,σRI=3,σSI=1,μTrue=100, and the standard deviation of the true values, σμ=1. The horizontal lines depict the five group means. Equation (8) implies that there is a need to partition error variance of d into “between” (B) and “within” (W) groups, as in(9)∑i=1g∑j=1ndij-d-2=∑i=1g∑j=1ndij-d-i2+n∑i=1gd-i-d-2=SSW+SSB.Figure 2
Values ofd=O-I in data simulated from (7) with n=10, g = 5, σRO=1,σSO=0.1,σRI=3,σSI=1,μTrue=100, and the SD of the true values σμ=1.


### 3.1. Grubbs’ Estimator as an Example of Top-Down UQ to EstimateσRO2,σSO2,σRI2,σSI2

Standard random-effects ANOVA [4] requires repeated measurements on some items in order to estimate σRI2 and then σSI2 in data assumed to be produced by a model such as (8). However, for most (O,I) data, repeated measurements of the same item are not available, so this section describes Grubbs’ estimator. Grubbs’ estimator was developed for situations in which more than one measurement method is applied to each of multiple test items (which may contain different material amounts), but there is no replication of measurements by any of the methods. Grubbs’ estimator will be described for additive measurement error models; a new version of Grubbs’ estimator for multiplicative error models is described in [1]. Note that the variance σRI2 of the random error variance component RIij includes the effects of “item-specific” bias (see Section 2, [7, 28]), which could not be estimated if available data were only from repeated measurements of the same or very similar items. Note also that Grubbs’ estimator does not consider the possibility of falsification by the operator, so it is intended to be applied to paired (O,I) data that has no evidence of falsification.The basis of Grubbs’ estimator within one group to estimateσRO2 and σRI2 is that the covariance between operator and inspector measurements equals the variance of the true item masses, σμ2, while the variance of I, σI2, conditional on the value of S is given by σI2=σμ2+σRI2. Therefore, the sample covariance within a single inspection period between operator and inspector measurements can be subtracted from the sample variance of the inspector measurements to estimate σRI2 (and similarly for estimating σRO2). That is, within one inspection period (group) (lowercase i(o) denotes the observed values of I(O)), Grubbs’ estimator is given by(10)σ^RI2=1n-1∑j=1nij-i-2-∑j=1noj-o-ij-i-.The estimates from (10) from each of the g groups are averaged to get the final estimate of the inspector’s random error variance, and similarly, the estimate of σμ2 is the average of the sample covariances σ^μ2=∑j=1noj-o-ij-i-/(n-1) computed within each group.To estimateσSI2 in (7), a minor extension of standard random-effects ANOVA to account for σμ2 shows that E{∑j=1gn(I-j-I--)2/(g-1)}=σRI2+σμ2+nσSI2, so a method-of-moments-based estimate of σSI2 is σ^SI2=(∑j=1g(i-j-i--)2)/(g-1)-(σ^RI2+σ^μ2)/n. There is no guarantee that σ^RI2 or σ^SI2 are nonnegative, but the corresponding true quantities are nonnegative (i.e., σRI2≥0 and σSI2≥0), so constrained versions of Grubbs-based and ANOVA-based estimators can be used (Section 3.2, [7, 30]).Grubbs showed (and simulation inR [25] also verified) that his estimator for σRI2 has variance σσ^RI22=2σRI4/(n-1)+(1/n-1)(σRO2σRI2+σRO2σμ2+σRI2σμ2), which is relatively large, particularly if σμ2 is comparable in magnitude to σRI2 [7] and/or n is small, so [13] proposed an option to mitigate the impact of σμ2 on σσ^RI22. The method in [13] relied on knowing the value, or approximate value of σRI2/σμ2, and studied the sensitivity to misspecifying the ratio σRI2/σμ2. The Bayesian option in Section 3.2 specifies a probability distribution for σRI2/σμ2 prior to observing the (O,I) data, as an example of enforcing nonnegativity constraints and including prior information, such as information from bottom-up UQ regarding σμ2,σRI2, and/or σSI2, and similarly for the operator or for variance ratios such as σRI2/σRO2.A measurement error model that can be used in top-down UQ for the type of data in Figure2 was given in (8). In the case of inverse regression as a type of bottom-up UQ, using (5), one can modify the top-down error model of (8) to(11)Iij=μij+S1Ii+S2Ii+RIij,where R~N0,σR2 (σR2 is random error variance that includes the effects of errors in X and Y); S1=α^0-α0 (an additive subcomponent of systematic error; see below), and S2=(α^1-α1)(YTest-y-) (a multiplicative subcomponent of systematic error).Assume that inverse regression is performed using mean-centered data (x~=x-x-train and y~=y-y-train), so cov(α^0,α^1)=0, which simplifies interpretation. The equation S1=α^0-α0 is then interpreted to mean that S1 has mean zero and variance σR2/n. Note that one cannot use paired (O,I) data to estimate σS12 and σS22, because the effects of σS12 and σS22 are confounded. However, simulation such as that used to produce Figure 1 provides a way to perform bottom-up prediction of top-down SD (or RSD) that can be compared to the SD (or RSD) observed in top-down evaluations. If the bottom-up predicted SD agrees well with that observed in top-down (O,I) data, then bottom-up UQ via simulation (because analytical approximations such as those plotted in Figure 1 in the calibration context have been shown to not be sufficiently accurate) can separately estimate σS12 and σS22, or, if fitting a zero-intercept model, (8) could be modified to a multiplicative error model, Iij=μij(1+SIi+RIij), where SI=α^1-α1YTest-y-, again with the variance of SI estimated by simulation [7, 17, 31].


### 3.2. New Bayesian Approach to Grubbs-Type Estimation

Recall that the variance of Grubbs’ estimator can be large, so [7, 31] review alternatives to Grubbs’ estimator based on constrained optimization, such as Jaech’s [31] constrained maximum likelihood estimator (CMLE), which assumes that the random and systematic measurement errors are normally distributed. Also, although the impact of σμ2 can be relatively large on Grubbs’ estimator, versions of Grubbs’ estimators that are constrained so that σ^RI2+σ^RO2=σ^d2, where σ^d2 is the sample variance of the differences d=O-I, have exhibited lower RMSE than Grubbs’ estimator in limited testing to date. Any estimator of σRO2,σRI2,σSO2, and σSI2 must be accompanied by its respective uncertainty. The uncertainty in CMLE or constrained least squares estimators [32] can be approximated using approximate analytical results and asymptotic results or by resampling methods such as the bootstrap. The quality of such approximations is not yet known; so a Bayesian alternative is presented here which does not rely on such approximations or the bootstrap for assessing uncertainty in σ^RO2, σ^RI2, σ^SO2, or σ^SI2.Another reason to consider the Bayesian alternative is that Grubbs’ estimator exhibits dependence on the underlying measurement error distribution. Figure3 plots the estimated probability density for Grubbs’ estimator for σ^RO, when the underlying random error distribution is either the normal, lognormal, or generalized lambda with thin or thick tails, for the relatively large sample size of 5 groups and 10 measurements per group. The probability density for σ^RO was estimated using a kernel-density estimator in R [25]. There is a relatively large uncertainty in σ^RO (Section 4 evaluates one impact of uncertainty in estimated RSDs), with an RSD in σ^RO of 11%, 8%, 13%, and 37%, respectively, for the four distributions in Figure 3, and an RSD in σ^SO of approximately 50% for all four distributions; also, for 2 groups and 5 measurements per group, the RSD in σ^RO is approximately 50% for all four distributions. One implication of Figure 3 is that uncertainties in σ^RO2,σ^RI2,σ^SO2, or σ^SI2 depend on the error distributions. An advantage of the Bayesian option is that the width of the Bayesian posterior adjusts to accommodate nonnormal underlying distributions [7].Figure 3
The estimated probability density of Grubbs’ estimator forσRO when the random errors have either a normal, lognormal, or generalized lambda (thin tail or thick tail) distribution. The true value of σRO is 0.0111, indicated by the vertical line.One option to improve Grubbs’ estimator is to impose constraints, such asσRO≤σRI, or, more flexibly, by assigning a prior probability to the ratio σRO/σRI which puts most of the probability on values less than 1. The uncertainty in constrained estimators is simple to estimate in a Bayesian approach and is often difficult to estimate in a non-Bayesian framework. Within an inspection group for fixed SO and SI, the paired (O,I) data from (8) has a bivariate distribution with covariance matrix ΣW=σμ2+σRO2σμ2σμ2σμ2+σRI2 and [33] provided a Bayesian approach to estimate σμ2, σRO2, and σRI2, assuming a bivariate normal likelihood, without imposing constraints on any of the variance ratios. In principle, one could extend the Bayesian approach in [33] to allow for a nonnormal likelihood and/or to allow for constraints on any of the variance ratios. In practice, such extensions are technically difficult and rarely attempted.In any Bayesian approach, prior information regarding the sizes or relative sizes ofσμ2, σRO2, and σRI2 must be provided. If the prior is “conjugate” for the likelihood, then the posterior is in the same likelihood family as the prior, in which case analytical methods are available to compute posterior prediction intervals for quantities of interest, so that a wide variety of priors and likelihoods can be accommodated; modern Bayesian methods do not rely on conjugate priors but use numerical methods to obtain samples of σμ2, σRO2, and σRI2 from their approximate posterior distributions [34]. For numerical methods such as Markov Chain Monte Carlo, the user must specify a prior distribution for σμ2, σRO2, and σRI2 and a likelihood (which need not be normal). The Bayesian approach presented next is approximate Bayesian computation (ABC), which does not require a known likelihood for the data and can accommodate constraints on variances and/or ratios of variances by choice of the prior distributions.The “output” of any Bayesian analysis is the posterior distribution and so the output of ABC is an estimate of the posterior distributions ofσμ2, σRO2, and σRI2. No matter what type of Bayesian approach is used, a well-calibrated Bayesian approach satisfies several requirements. The requirement of interest here is that, in repeated applications of ABC, approximately 95% of the middle 95% of the posterior distribution for each of σμ2, σRO2, and σRI2 should contain the respective true values.In ABC, one assumes that a model has input parametersθ and outputs data yM=y(θ) (M for “model”) and that there is corresponding observed real data yobs. Here, the model is (8), which specifies how to generate synthetic O and I data and does require a likelihood; however, the true likelihood used to generate the data need not be known to the user. Synthetic data is generated from the model for many trial values of θ, and trial θ values are accepted as contributing to the estimated posterior distribution for θyobs if the distance D(yobs,y(θ)) between yobs and yM=y(θ) is reasonably small. Alternatively, for most applications, it is necessary to reduce the dimension of yobs to a small set of summary statistics S and instead accept trial values of θ if D(S(yobs), S(y(θ)))<T, where T is a user-chosen threshold. Here, yobs is the paired (O,I) data in each inspection group, and S(yobs) includes within- and between-groups sums of squares and within-group covariance between O and I. Specifically, recall that the estimator of σRI2 in (8) is σ^RI2=1/n-1∑j=1nij-i-2-∑j=1noj-o-ij-i- and that a method-of-moments-based estimate of σSI2 is σ^SI2=(∑j=1g(i-j-i--)2)/(g-1)-(σ^RI2+σ^μ2)/n. And σμ2 can be estimated using σ^μ2=∑j=1n(oj-o-)ij-i-. The quantities σ^RI2,σ^SI2,σ^RO2,σ^SO2, and σ^μ2 are therefore good choices for summary statistics to be used for ABC. Because trial values of θ are accepted if D(S(yobs), S(y(θ)))<T, an approximation error to the posterior distribution arises which several ABC options attempt to mitigate. Such options involve weighting the accepted θ values by the actual distance D(S(yobs), S(y(θ))) (abctools in R [25]). As an aside, if the error model is multiplicative rather than additive, expressions given in [1] can be used as summary statistics.To summarize, ABC consists of three steps: (1) sample parameter values from their prior distribution pprior(θ); (2) for each simulated value of θ in (1), simulate data from (8); (3) accept a fraction of the sampled prior values in (1) by checking whether the summary statistics computed from the data in (2) satisfy D(S(yobs), S(y(θ)))<T. If desired, aiming to improve the approximation to the posterior, adjust the accepted θ values on the basis of the actual D(yobs,θ) value. ABC requires the user to make three choices: the summary statistics, the threshold T, and the measure of distance d. Reference [35] introduced a method to choose summary statistics that use the estimated posterior means of the parameters based on pilot simulation runs. Reference [36] used an estimate of the change in posterior pposterior(θ) when candidate summary statistics are added to the current set of summary statistics. Reference [37] illustrated a method to evaluate whether a candidate set of summary statistics leads to a well-calibrated posterior in the same sense that is used in this paper; that is, nominal posterior probability intervals should have approximately the same actual coverage probability.To illustrate application of ABC to top-down UQ, simulations using (8) were performed. Recall that an additive model is a reasonable approximation to a multiplicative model if the error variances are small and the data is analyzed on a log scale or if there are effects in addition to calibration which impact the random and systematic errors, such as random changes in background count rates which are not adjusted for. Simulations from a multiplicative error model were also investigated, with good results, such as those described next for the additive model; the summary statistics for ABC to use Grubbs-type estimators for a multiplicative model are given in [1, 7].The simulations were performed inR using three steps. In the first step, ABC requires a training collection of parameter values and corresponding summary statistics for each of many simulations. So, in each of 105 simulations, the values for σRO,σSO,σRI,σSI,σμ were sampled from their respective prior probability densities. In the second step, (8) was used to generate I and O data. In the third step, the expressions for σ^RO2,σ^SO2,σ^RI2,σ^SI2, and σ^μ2 given above were used as summary statistics, resulting in a parameter matrix and corresponding summary statistics matrix, each having 105 rows and 5 columns. Then, in a separate set of 105 simulations, the same first and second steps were repeated, and for each simulated set of parameters and summary statistics, the parameter and summary statistics matrices from the third step in training were used in the abc function in the abctools package, using an acceptance fraction of 0.01 (meaning that 1% of the trial values for the true parameters were accepted), to produce an approximate posterior for each of the five parameters. This posterior can be used to assess the actual coverage probability, for example, of an interval that contains 95% of the posterior probability.It was found [7] that the actual coverages for σRO,σSO,σRI,σSI were essentially the same (to within simulation uncertainty) as the nominal coverages, at 90%, 95%, and 99% probabilities, for a normal distribution and all of the nonnormal distributions investigated (uniform, gamma, lognormal, beta,t, and generalized lambda with thick or thin tails), each having the same σRO,σSO,σRI,σSI values. In addition, when a normal distribution was used to train ABC and any of the evaluated nonnormal distributions were used to test ABC, very nearly the same actual coverages (to within approximately ±0.01) were obtained. As another check of robustness, one prior distribution was used for training ABC and a prior that was wider by a factor of 1.3 was used in testing. In that case, the actual coverages dropped from approximately 0.95 to approximately 0.90, so this implementation is less robust to using a wider prior in testing than in training. Also, the RMSE of the ABC estimator was compared to each of several non-Bayesian constrained estimators that are currently being evaluated. Not surprisingly, provided that the prior used in training was approximately the same as that used in testing, the ABC estimator had lowest RMSE. The intent here is not to make any general RMSE performance claims; instead, these results provide a second indication that the ABC implementation is well calibrated, in the sense that if the assumed prior equals the true prior then the RMSE of the ABC estimator is highly competitive with that of the non-Bayesian estimator and in the sense that the width of the posterior accurately describes uncertainty. Finally, a well-calibrated Bayesian analysis allows one to evaluate the effect of increasing sample size on uncertainty in the estimated values of σRO,σSO,σRI,σSI.To illustrate ABC output, Figure4 plots the prior density and the estimated posterior density for σRI and σSI for g = 3 and n=3 and for g = 10 and n=10. Because the posterior densities are well calibrated, they can be used to reliably assess whether top-down estimates of σRI and σSI are in agreement to within their respective uncertainties of the corresponding bottom-up estimates of σRI and σSI from Section 2. The priors used in Figure 4 are wide, with no relative information regarding variance ratios assumed, and the parameter σμ is assigned a prior with a mean value of 0.10, with a relative standard deviation in σμ of 10%. And Figure 4 shows that, in this case, having only g = 3 and n=3 per group does not lead to a narrow posterior, but, with g = 10 and n=10, the posterior is fairly narrow.Figure 4
Prior and posterior probability densities for data from realization of (8) for σRI and σSI for g = 3 and n=3 and for g = 10 and n=10. The true values are σRI = 2.3 (% relative) and σSI = 1.2 (% relative).
(a)
σ
R
(b)
σ
S


## 3.1. Grubbs’ Estimator as an Example of Top-Down UQ to EstimateσRO2,σSO2,σRI2,σSI2

Standard random-effects ANOVA [4] requires repeated measurements on some items in order to estimate σRI2 and then σSI2 in data assumed to be produced by a model such as (8). However, for most (O,I) data, repeated measurements of the same item are not available, so this section describes Grubbs’ estimator. Grubbs’ estimator was developed for situations in which more than one measurement method is applied to each of multiple test items (which may contain different material amounts), but there is no replication of measurements by any of the methods. Grubbs’ estimator will be described for additive measurement error models; a new version of Grubbs’ estimator for multiplicative error models is described in [1]. Note that the variance σRI2 of the random error variance component RIij includes the effects of “item-specific” bias (see Section 2, [7, 28]), which could not be estimated if available data were only from repeated measurements of the same or very similar items. Note also that Grubbs’ estimator does not consider the possibility of falsification by the operator, so it is intended to be applied to paired (O,I) data that has no evidence of falsification.The basis of Grubbs’ estimator within one group to estimateσRO2 and σRI2 is that the covariance between operator and inspector measurements equals the variance of the true item masses, σμ2, while the variance of I, σI2, conditional on the value of S is given by σI2=σμ2+σRI2. Therefore, the sample covariance within a single inspection period between operator and inspector measurements can be subtracted from the sample variance of the inspector measurements to estimate σRI2 (and similarly for estimating σRO2). That is, within one inspection period (group) (lowercase i(o) denotes the observed values of I(O)), Grubbs’ estimator is given by(10)σ^RI2=1n-1∑j=1nij-i-2-∑j=1noj-o-ij-i-.The estimates from (10) from each of the g groups are averaged to get the final estimate of the inspector’s random error variance, and similarly, the estimate of σμ2 is the average of the sample covariances σ^μ2=∑j=1noj-o-ij-i-/(n-1) computed within each group.To estimateσSI2 in (7), a minor extension of standard random-effects ANOVA to account for σμ2 shows that E{∑j=1gn(I-j-I--)2/(g-1)}=σRI2+σμ2+nσSI2, so a method-of-moments-based estimate of σSI2 is σ^SI2=(∑j=1g(i-j-i--)2)/(g-1)-(σ^RI2+σ^μ2)/n. There is no guarantee that σ^RI2 or σ^SI2 are nonnegative, but the corresponding true quantities are nonnegative (i.e., σRI2≥0 and σSI2≥0), so constrained versions of Grubbs-based and ANOVA-based estimators can be used (Section 3.2, [7, 30]).Grubbs showed (and simulation inR [25] also verified) that his estimator for σRI2 has variance σσ^RI22=2σRI4/(n-1)+(1/n-1)(σRO2σRI2+σRO2σμ2+σRI2σμ2), which is relatively large, particularly if σμ2 is comparable in magnitude to σRI2 [7] and/or n is small, so [13] proposed an option to mitigate the impact of σμ2 on σσ^RI22. The method in [13] relied on knowing the value, or approximate value of σRI2/σμ2, and studied the sensitivity to misspecifying the ratio σRI2/σμ2. The Bayesian option in Section 3.2 specifies a probability distribution for σRI2/σμ2 prior to observing the (O,I) data, as an example of enforcing nonnegativity constraints and including prior information, such as information from bottom-up UQ regarding σμ2,σRI2, and/or σSI2, and similarly for the operator or for variance ratios such as σRI2/σRO2.A measurement error model that can be used in top-down UQ for the type of data in Figure2 was given in (8). In the case of inverse regression as a type of bottom-up UQ, using (5), one can modify the top-down error model of (8) to(11)Iij=μij+S1Ii+S2Ii+RIij,where R~N0,σR2 (σR2 is random error variance that includes the effects of errors in X and Y); S1=α^0-α0 (an additive subcomponent of systematic error; see below), and S2=(α^1-α1)(YTest-y-) (a multiplicative subcomponent of systematic error).Assume that inverse regression is performed using mean-centered data (x~=x-x-train and y~=y-y-train), so cov(α^0,α^1)=0, which simplifies interpretation. The equation S1=α^0-α0 is then interpreted to mean that S1 has mean zero and variance σR2/n. Note that one cannot use paired (O,I) data to estimate σS12 and σS22, because the effects of σS12 and σS22 are confounded. However, simulation such as that used to produce Figure 1 provides a way to perform bottom-up prediction of top-down SD (or RSD) that can be compared to the SD (or RSD) observed in top-down evaluations. If the bottom-up predicted SD agrees well with that observed in top-down (O,I) data, then bottom-up UQ via simulation (because analytical approximations such as those plotted in Figure 1 in the calibration context have been shown to not be sufficiently accurate) can separately estimate σS12 and σS22, or, if fitting a zero-intercept model, (8) could be modified to a multiplicative error model, Iij=μij(1+SIi+RIij), where SI=α^1-α1YTest-y-, again with the variance of SI estimated by simulation [7, 17, 31].


## 3.2. New Bayesian Approach to Grubbs-Type Estimation

Recall that the variance of Grubbs’ estimator can be large, so [7, 31] review alternatives to Grubbs’ estimator based on constrained optimization, such as Jaech’s [31] constrained maximum likelihood estimator (CMLE), which assumes that the random and systematic measurement errors are normally distributed. Also, although the impact of σμ2 can be relatively large on Grubbs’ estimator, versions of Grubbs’ estimators that are constrained so that σ^RI2+σ^RO2=σ^d2, where σ^d2 is the sample variance of the differences d=O-I, have exhibited lower RMSE than Grubbs’ estimator in limited testing to date. Any estimator of σRO2,σRI2,σSO2, and σSI2 must be accompanied by its respective uncertainty. The uncertainty in CMLE or constrained least squares estimators [32] can be approximated using approximate analytical results and asymptotic results or by resampling methods such as the bootstrap. The quality of such approximations is not yet known; so a Bayesian alternative is presented here which does not rely on such approximations or the bootstrap for assessing uncertainty in σ^RO2, σ^RI2, σ^SO2, or σ^SI2.Another reason to consider the Bayesian alternative is that Grubbs’ estimator exhibits dependence on the underlying measurement error distribution. Figure3 plots the estimated probability density for Grubbs’ estimator for σ^RO, when the underlying random error distribution is either the normal, lognormal, or generalized lambda with thin or thick tails, for the relatively large sample size of 5 groups and 10 measurements per group. The probability density for σ^RO was estimated using a kernel-density estimator in R [25]. There is a relatively large uncertainty in σ^RO (Section 4 evaluates one impact of uncertainty in estimated RSDs), with an RSD in σ^RO of 11%, 8%, 13%, and 37%, respectively, for the four distributions in Figure 3, and an RSD in σ^SO of approximately 50% for all four distributions; also, for 2 groups and 5 measurements per group, the RSD in σ^RO is approximately 50% for all four distributions. One implication of Figure 3 is that uncertainties in σ^RO2,σ^RI2,σ^SO2, or σ^SI2 depend on the error distributions. An advantage of the Bayesian option is that the width of the Bayesian posterior adjusts to accommodate nonnormal underlying distributions [7].Figure 3
The estimated probability density of Grubbs’ estimator forσRO when the random errors have either a normal, lognormal, or generalized lambda (thin tail or thick tail) distribution. The true value of σRO is 0.0111, indicated by the vertical line.One option to improve Grubbs’ estimator is to impose constraints, such asσRO≤σRI, or, more flexibly, by assigning a prior probability to the ratio σRO/σRI which puts most of the probability on values less than 1. The uncertainty in constrained estimators is simple to estimate in a Bayesian approach and is often difficult to estimate in a non-Bayesian framework. Within an inspection group for fixed SO and SI, the paired (O,I) data from (8) has a bivariate distribution with covariance matrix ΣW=σμ2+σRO2σμ2σμ2σμ2+σRI2 and [33] provided a Bayesian approach to estimate σμ2, σRO2, and σRI2, assuming a bivariate normal likelihood, without imposing constraints on any of the variance ratios. In principle, one could extend the Bayesian approach in [33] to allow for a nonnormal likelihood and/or to allow for constraints on any of the variance ratios. In practice, such extensions are technically difficult and rarely attempted.In any Bayesian approach, prior information regarding the sizes or relative sizes ofσμ2, σRO2, and σRI2 must be provided. If the prior is “conjugate” for the likelihood, then the posterior is in the same likelihood family as the prior, in which case analytical methods are available to compute posterior prediction intervals for quantities of interest, so that a wide variety of priors and likelihoods can be accommodated; modern Bayesian methods do not rely on conjugate priors but use numerical methods to obtain samples of σμ2, σRO2, and σRI2 from their approximate posterior distributions [34]. For numerical methods such as Markov Chain Monte Carlo, the user must specify a prior distribution for σμ2, σRO2, and σRI2 and a likelihood (which need not be normal). The Bayesian approach presented next is approximate Bayesian computation (ABC), which does not require a known likelihood for the data and can accommodate constraints on variances and/or ratios of variances by choice of the prior distributions.The “output” of any Bayesian analysis is the posterior distribution and so the output of ABC is an estimate of the posterior distributions ofσμ2, σRO2, and σRI2. No matter what type of Bayesian approach is used, a well-calibrated Bayesian approach satisfies several requirements. The requirement of interest here is that, in repeated applications of ABC, approximately 95% of the middle 95% of the posterior distribution for each of σμ2, σRO2, and σRI2 should contain the respective true values.In ABC, one assumes that a model has input parametersθ and outputs data yM=y(θ) (M for “model”) and that there is corresponding observed real data yobs. Here, the model is (8), which specifies how to generate synthetic O and I data and does require a likelihood; however, the true likelihood used to generate the data need not be known to the user. Synthetic data is generated from the model for many trial values of θ, and trial θ values are accepted as contributing to the estimated posterior distribution for θyobs if the distance D(yobs,y(θ)) between yobs and yM=y(θ) is reasonably small. Alternatively, for most applications, it is necessary to reduce the dimension of yobs to a small set of summary statistics S and instead accept trial values of θ if D(S(yobs), S(y(θ)))<T, where T is a user-chosen threshold. Here, yobs is the paired (O,I) data in each inspection group, and S(yobs) includes within- and between-groups sums of squares and within-group covariance between O and I. Specifically, recall that the estimator of σRI2 in (8) is σ^RI2=1/n-1∑j=1nij-i-2-∑j=1noj-o-ij-i- and that a method-of-moments-based estimate of σSI2 is σ^SI2=(∑j=1g(i-j-i--)2)/(g-1)-(σ^RI2+σ^μ2)/n. And σμ2 can be estimated using σ^μ2=∑j=1n(oj-o-)ij-i-. The quantities σ^RI2,σ^SI2,σ^RO2,σ^SO2, and σ^μ2 are therefore good choices for summary statistics to be used for ABC. Because trial values of θ are accepted if D(S(yobs), S(y(θ)))<T, an approximation error to the posterior distribution arises which several ABC options attempt to mitigate. Such options involve weighting the accepted θ values by the actual distance D(S(yobs), S(y(θ))) (abctools in R [25]). As an aside, if the error model is multiplicative rather than additive, expressions given in [1] can be used as summary statistics.To summarize, ABC consists of three steps: (1) sample parameter values from their prior distribution pprior(θ); (2) for each simulated value of θ in (1), simulate data from (8); (3) accept a fraction of the sampled prior values in (1) by checking whether the summary statistics computed from the data in (2) satisfy D(S(yobs), S(y(θ)))<T. If desired, aiming to improve the approximation to the posterior, adjust the accepted θ values on the basis of the actual D(yobs,θ) value. ABC requires the user to make three choices: the summary statistics, the threshold T, and the measure of distance d. Reference [35] introduced a method to choose summary statistics that use the estimated posterior means of the parameters based on pilot simulation runs. Reference [36] used an estimate of the change in posterior pposterior(θ) when candidate summary statistics are added to the current set of summary statistics. Reference [37] illustrated a method to evaluate whether a candidate set of summary statistics leads to a well-calibrated posterior in the same sense that is used in this paper; that is, nominal posterior probability intervals should have approximately the same actual coverage probability.To illustrate application of ABC to top-down UQ, simulations using (8) were performed. Recall that an additive model is a reasonable approximation to a multiplicative model if the error variances are small and the data is analyzed on a log scale or if there are effects in addition to calibration which impact the random and systematic errors, such as random changes in background count rates which are not adjusted for. Simulations from a multiplicative error model were also investigated, with good results, such as those described next for the additive model; the summary statistics for ABC to use Grubbs-type estimators for a multiplicative model are given in [1, 7].The simulations were performed inR using three steps. In the first step, ABC requires a training collection of parameter values and corresponding summary statistics for each of many simulations. So, in each of 105 simulations, the values for σRO,σSO,σRI,σSI,σμ were sampled from their respective prior probability densities. In the second step, (8) was used to generate I and O data. In the third step, the expressions for σ^RO2,σ^SO2,σ^RI2,σ^SI2, and σ^μ2 given above were used as summary statistics, resulting in a parameter matrix and corresponding summary statistics matrix, each having 105 rows and 5 columns. Then, in a separate set of 105 simulations, the same first and second steps were repeated, and for each simulated set of parameters and summary statistics, the parameter and summary statistics matrices from the third step in training were used in the abc function in the abctools package, using an acceptance fraction of 0.01 (meaning that 1% of the trial values for the true parameters were accepted), to produce an approximate posterior for each of the five parameters. This posterior can be used to assess the actual coverage probability, for example, of an interval that contains 95% of the posterior probability.It was found [7] that the actual coverages for σRO,σSO,σRI,σSI were essentially the same (to within simulation uncertainty) as the nominal coverages, at 90%, 95%, and 99% probabilities, for a normal distribution and all of the nonnormal distributions investigated (uniform, gamma, lognormal, beta,t, and generalized lambda with thick or thin tails), each having the same σRO,σSO,σRI,σSI values. In addition, when a normal distribution was used to train ABC and any of the evaluated nonnormal distributions were used to test ABC, very nearly the same actual coverages (to within approximately ±0.01) were obtained. As another check of robustness, one prior distribution was used for training ABC and a prior that was wider by a factor of 1.3 was used in testing. In that case, the actual coverages dropped from approximately 0.95 to approximately 0.90, so this implementation is less robust to using a wider prior in testing than in training. Also, the RMSE of the ABC estimator was compared to each of several non-Bayesian constrained estimators that are currently being evaluated. Not surprisingly, provided that the prior used in training was approximately the same as that used in testing, the ABC estimator had lowest RMSE. The intent here is not to make any general RMSE performance claims; instead, these results provide a second indication that the ABC implementation is well calibrated, in the sense that if the assumed prior equals the true prior then the RMSE of the ABC estimator is highly competitive with that of the non-Bayesian estimator and in the sense that the width of the posterior accurately describes uncertainty. Finally, a well-calibrated Bayesian analysis allows one to evaluate the effect of increasing sample size on uncertainty in the estimated values of σRO,σSO,σRI,σSI.To illustrate ABC output, Figure4 plots the prior density and the estimated posterior density for σRI and σSI for g = 3 and n=3 and for g = 10 and n=10. Because the posterior densities are well calibrated, they can be used to reliably assess whether top-down estimates of σRI and σSI are in agreement to within their respective uncertainties of the corresponding bottom-up estimates of σRI and σSI from Section 2. The priors used in Figure 4 are wide, with no relative information regarding variance ratios assumed, and the parameter σμ is assigned a prior with a mean value of 0.10, with a relative standard deviation in σμ of 10%. And Figure 4 shows that, in this case, having only g = 3 and n=3 per group does not lead to a narrow posterior, but, with g = 10 and n=10, the posterior is fairly narrow.Figure 4
Prior and posterior probability densities for data from realization of (8) for σRI and σSI for g = 3 and n=3 and for g = 10 and n=10. The true values are σRI = 2.3 (% relative) and σSI = 1.2 (% relative).
(a)
σ
R
(b)
σ
S


## 4. Application of RSD Estimates: Statistical Testing of Materials Accounting Data

### 4.1. Sequential Statistical Testing of MB Sequences

Recall that NMA evaluates one or more MBs, where the MB is defined for balance periodj as MBj=(Ij-1+Tin,j-Tout,j)-Ij [38, 39]. Typically, many measurements are combined to estimate the terms Tin,Ibegin,Tout, and Iend in the MB; therefore, the central limit effect and years of experience suggest that MBs in most facilities will be approximately normally distributed with the mean equal to the true NM loss λj  and standard deviation σj, which is expressed as Xj~N(λj,σj), where X denotes the MB and the notation σj is a shortened version of σMB,j. A sequence of n MBs will be assumed to have approximately a multivariate normal distribution [38–43],  X1,X2,…,Xn~MVN(λ,Σ), where then-by-n covariance matrix is(12)Σ=σ12σ122⋯σ1n2σ212σ22⋯σ2n2⋮σn12σn22⋯σn2.The magnitude ofσj determines the amount of NM loss, λ=∑k=1nλk, which leads to high detection probability (DP). For example, suppose the facility tests for NM loss only, not for NM gain, and assume that Xj~N(λj,σj) is an adequate model. Then, if a false alarm probability (FAP) of α = 0.05 is desired, the alarm threshold is 1.65σj. In the case of testing for loss only, it follows that the loss detection probability 1-β for λ = 3.3σ and 1-β > 0.95 if λ> 3.3σj, where β is the nondetection (false negative) probability. The factor 3.3 arises from symmetry of the Gaussian distribution, requiring α = β = 0.05, and the fact that 1.65 = 3.3/2 is the 0.95 quantile of the N(0,1) distribution. One common goal is for DP = 1-β to be at least 0.95 if λ ≥ 1 SQ (significant quantity, which is 8 kg for Pu), which is accomplished if σj≤ SQ/3.3. If σj > SQ/3.3, this can be mitigated by reducing measurement errors to achieve σj≤ SQ/3.3 (if feasible) and/or by closing the balances more frequently, so there is less nuclear material transferred per balance period, which reduces σj [38, 39]. The DP of other safeguards measures such as enhanced containment and surveillance with smart cameras and/or remote radiation detection is difficult to quantify and is outside the scope of this paper. This section concludes with four remarks.Remark 1.
Large throughput facilities try to makeσj as small as reasonably possible and often try to keep σj small as a percent of throughput but cannot achieve σj≤ SQ/3.3. For example, suppose that there is a measurement error relative standard deviation of 0.5% of throughput. And suppose that the FAP/DP goals are α = 0.05 and 1-β = 0.95 and annual throughput is 100 SQ. Then, σj = 0.5 SQ = SQ/2 > SQ/3.3, so protracted diversion of 1 SQ over one year will not have a high DP. Therefore, one reason for frequent MB accounting is that abrupt diversion over hours or days is very likely to be detected [40]. As a complementary approach that is beyond the scope here, process monitoring [39] methods have the potential to detect off-normal facility operation that could misdirect NM to undeclared locations.Remark 2.
In the 1980s, some believed that a plant reporting a MB every 30 days would have higher DP than that same plant reporting a MB every year. However, [44] showed that, for optimal (from the diverter’s viewpoint) protracted diversion with the per-period loss being proportional to the row sums of the covariance matrix Σ of the MB sequence, annual MB testing has higher DP than monthly MB testing, and so, for such protracted diversion, “less frequent balance closure is better.” However, [44] conceded that NRTA has shorter detection times and higher DP against abrupt diversion. Publications [45–50] soon followed involving joint sequential tests: one tuned to detect protracted diversion and one more tuned for abrupt diversion. Such joint Page’s tests can be tuned to have high DP against abrupt loss while still having reasonably high DPs against protracted loss. Two types of combined Page’s tests are included in the simulation study in Section 4.3.Remark 3.
The assumptionX1,X2,…,Xn~MVN(λ,Σ) implies that Σ is known without estimation error. In practice, Σ is estimated using variance propagation applied to Xj=MBj=(Ij-1+Tin,j-Tout,j)-Ij and there will be estimation error in the estimate Σ^ [40].
The simulation study in Section4.3 includes a sensitivity study to assess the impact of estimation error Σ^ on the estimated false alarm probabilities (FAPs) and detection probabilities (DPs).Remark 4.
This section focuses on operator MB sequences. In international safeguards, the inspector randomly selects items to verify NM declarations made by the operator. The difference statistic,D=N∑j=1n((oj-ij)/n), defined as the average difference in the sample of sizen (extrapolated to the population of sizeN) between operator declarations (almost always based on operator measurements) and inspector measurements can be used as a test statistic, or the D statistic [2] can be used to compute the inspector MB statistic (or sequence). Inspector MB=MB-D, which could be analyzed using sequential statistical methods such as those in Section 4.3. Alternatively, one can test individually each of the n paired differences dj=oj-ij and the overall D statistic, and if none are found to be statistically significant, then the IAEA could rely on operator MB evaluation as in Section 4.3.


### 4.2. Propagation of Variance to EstimateΣ

EstimatingΣ is a key step required in frequent NMA. To illustrate, a simplified example model of a generic electrochemical facility with one input stream, one output stream, and one key inventory item will be used here [38]. Each measurement method is modeled here using a multiplicative measurement error model for the operator (O): Mi=μi(1+Si+Ri), with Si~N(0,δS2) and Ri~N(0,δR2), where Mi is the operator’s measured value of item i, μi is the true but unknown value of item i, Ri is a random error of item i, and Si is a short-term systematic error for item i. Then, the diagonal terms of Σ are calculated as(13)σi2=Tini2δin,R2+δin,S2+Touti2σout,R2+σout,S2+Ii2δinv,R2+Ii-12δinv,R2+Ii-Ii-12δinv,S2.The off-diagonal terms in Σ are calculated as(14)σij2=TiniTinjδin,S2+ToutiToutjδout,S2+IiIj+Ii-1Ij-1δinv,S2-IiIj-1δinv,S2+δinv,R2ifj-i=1-Ii-1Ijδinv,S2+δinv,R2ifi-j=1.In the last two terms, the random error of the inventory term is only applied if the condition is true. Reference [31] showed that the POV results for σi2 and σij2 are obtained by appropriate interpretation of GUM’s equation (1) in Section 2. For this simplified version of an example from [38], this leads to examples of 12-by-12 covariance matrices for monthly MBs over a one-year period; three example matrices Σ1, Σ2, and Σ3 are listed next.Three example covariance matrices,Σ1, Σ2, and Σ3, for a generic electrochemical facility, are given [38]:(15)1.00-0.480.010.010.01-0.481.00-0.480.010.010.01-0.481.00-0.480.010.010.01-0.481.00-0.480.010.010.01-0.481.00Equation (15) is Σ1, scaled to unit variance, only displaying 5 by 5 of the 12 by 12. This is the nominal case with 4 kg input per period and 40 kg inventory; δR = δS = 0.01.(16)1.000.170.330.330.330.171.000.170.330.330.330.171.000.170.330.330.330.171.000.170.330.330.330.171.00Equation (16) is Σ2, scaled to unit variance. This is the case with 4 kg input per period and 4 kg inventory.(17)1.000.580.670.670.670.581.000.580.670.670.670.581.000.580.670.670.670.581.000.580.670.670.670.581.00Equation (17) is Σ3, scaled to unit variance. This is the case with 4 kg input per period and 40 kg inventory and δS increased from 0.01 to 0.02.


### 4.3. Sequential Testing of Material Balance in Nuclear Material Accountancy

The assumptionX1,X2,…,Xn~MVN(λ,Σ) implies that Y=Σ-1/2X~MVN(Σ-1/2λ,I), where I is the identity matrix. The transform Y=Σ-1/2X is known in safeguards as the standardized independently transformed MUF (SITMUF, where MUF is another name for the MB), which is most conveniently computed using the Cholesky decomposition [43]. There are two main advantages of applying statistical tests to Y rather than to X. First, alarm thresholds depend only on the sequence length n forY and not on the form of the covariance matrix Σ. Because it is best to calculate thresholds using simulation, this is a logistic advantage. Second, the variance of the Y sequence decreases over time, so if a diversion occurs late in the analysis period, the DP is larger for the Y sequence than for the X sequence. Note that one cannot claim higher DP for the Y sequence than for the X sequence in general, because the true loss scenario is never known, and the DP can be larger for X than for Y for some loss scenarios, which is demonstrated in Section 4.The value ofYi can be calculated using Y=Σ-1/2X, but, more intuitively as the residual from the X sequence, Yj={Xj-E(Xj∣Xj-1,Xj-2,…,X1)}/σ~j, where the standard deviation σ~j is given by σ~j=σjj2-fΣ-1fT, where f=Σj,1:(j-1), the 1 to (j-1) entries in the jth row of Σ.Several reasonable statistical tests have been evaluated in [38, 41, 44–51] and are included in the simulation study in Section 4.4, including the following 13 tests:(1)
MUF test: this compares each MUF value to a threshold, which is the same as a Shewhart test in quality control (QC). The test alarms on periodj if MUFj/σj≥T for some threshold T(2)
SITMUF test: this compares each SITMUF valueY=Σ-1/2X to a threshold, which is the same as a Shewhart test in QC. The test alarms on period j if SITMUFj/σSITMUF,j≥T for some thresholdT(3)
Page’s test applied to MUF: Page’s test to test for loss is a sequence of sequential probability ratio tests, defined asPj=max⁡(0,Pj-1+xj/σj-k), where P0=0 [52]. The test alarms on period j if Pj>T. The parameter k is a control parameter that is optimal for detecting a shift from zero loss to loss λ if k=λ/2. The alarm thresholdT (usually denoted as h in literature on Page’s test, but this paper uses T for alarm threshold) is chosen so that the FAP per analysis period (usually one year) is 0.05 or whatever FAP is specified(4)
Page’s test applied to SITMUF: Page’s test to test for loss is a sequence of sequential probability ratio tests, asPj=max(0,Pj-1+yj-k). The alarm threshold T is chosen so that the FAP per analysis period (usually one year) is 0.05 or whatever FAP is specified(5)
Combined Page’s tests applied to MUF: the use of Page’s test with a large value ofk and small value of T has good DP for abrupt loss, and the use of Page’s test with a small value of k and large value of T has good DP for protracted loss. Therefore, a reasonable option is to use a combination of two Page’s tests, one with large k and one with smallk(6)
Applying combined Page’s tests to SITMUF(7)
CUMUF: at periodj, CUMUFj=∑i=1jxi is the sum of all MUF values from period 1 to j(8)
GEMUF: it has been shown that if the loss vectorλ is known, then the Neyman-Pearson test statistic is λTΣ-1x, which is known as a matched filter in some literature. The GEMUF statistic substitutes xT for λT, so GEMUF=xTΣ-1x. In simulation studies, λ is known, so the NP test statistic is useful for calculating the largest possible DP. The GE in GEMUF is a German abbreviation of Geschätzter, which means “estimated,” so GEMUF means estimated MUF, and GEMUF is the same as the Mahalanobis distance from the0 vector and Hotelling’s multivariate T statistic [51](9)
A nonsequential version of the Neyman-Pearson test,λTΣ-1x, is useful to calculate the largest possible DP for given Σ and λ. For completeness, four other combined tests are also considered(10)
SITMUF and CUMUF(11)
Page’s on SITMUF and CUMUF(12)
SITMUF, Page’s on SITMUF, and GEMUF(13)
Page’s on SITMUF, CUMUF, and GEMUFThis section concludes with four remarks.Remark 1 (SITMUF transform).
The SITMUF transform is recommended for two reasons. First, simulation is typically used to select alarm thresholds, and it is convenient to always work on the same scale when selecting alarm thresholds, so the fact thatY=Σ-1/2X~MVN(Σ-1/2λ,I) is convenient. Note that alarm thresholds could be selected on the basis of exact or approximate analytical results for some, but not all, of the tests. For example, there are approximate expressions for T andk [53]. Second, the standard deviation σ~j is given by σ~j=σj2-fΣ-1fT, where f=Σj,1:(j-1), the 1 to (j-1) entries in thejth row of Σ, so the standard deviation of the MUF residuals decreases in the later periods. Therefore, the independence transform is analogous to a bias adjustment, leading to smaller prediction variance in later periods, which tends to increase the DP for SITMUF compared to MUF (there are exceptions where the DP for MUF is larger than the DP for SITMUF; see Section 4.4, DP results).Remark 2 (choosing thresholds).
Thresholds can be chosen in many ways and can be assumed to be constant for each period or not. Therefore, simulation DP results in Section4.4 are not claimed to be optimal but are example DP results.Remark 3 (performance criteria).
The main performance criterion for comparing tests is the DP. But the average time to detection and robustness to misspecifying the covariance matrixΣ are also important.Remark 4.
There are other tests in the literature, including the power one and scan tests [38, 41–43].


### 4.4. Simulation Study

Example DP results are in Tables1–3 for Σ1, Σ2, and Σ3, respectively, using 105 simulations (so are repeatable to within ±0.01) to choose alarm thresholds and to estimate DPs. Tables 1–3 indicate that, as expected, there is no overall best test. However, Page’s test on SITMUF or a combined Page’s test on SITMUF is often among the best performers. Whether Page’s test on SITMUF has larger DP than Page’s test on MUF depends on the covariance matrix.Table 1
DPs forΣ1. Boldface entries have the largest DPs (excluding NP) for the respective column.
DP
Loss 1(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)/10
Loss 2(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3)
Loss 3(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)
Loss 4(0,1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
MUF
0.07
0.58, 0.58, 0.58
0.15
0.19
SITMUF
0.20
0.65, 0.89, 0.81
0.82
0.58
Page on MUF
0.18
0.75, 0.85, 0.74
0.76
0.58
Page on SITMUF
0.71
0.82, 0.99, 0.56
1.0
0.99
Combined Page MUF
0.74
0.83, 0.99, 0.80
1.0
0.99
Combined Page, SITMUF
0.74
0.83, 0.98, 0.80
1.0
0.99
CUMUF
0.26
0.91, 0.64, 0.18
0.50
0.65
GEMUF
0.13
0.87, 0.70, 0.11
0.57
0.48
NP
0.82
0.99, 1.0, 0.99
1.0
1.0
SITMUF and CUMUF
0.23
0.80, 0.87, 0.80
0.81
0.62
SITMUF, Page, CUMUF
0.28
0.85, 0.85, 0.75
0.76
0.58
SITMUF, Page, GEMUF
0.20
0.84, 0.84, 0.72
0.74
0.61
Page, CUMUF, GEMUF
0.23
0.85, 0.75, 0.75
0.76
0.58Table 2
DPs forΣ2.
DP
Loss 1(0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3)
Loss 2(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3)
Loss 3(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)
Loss 4(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
MUF
0.07
0.60, 0.60, 0.60
0.18
0.18
SITMUF
0.06
0.62, 0.76, 0.80
0.14
0.14
Page on MUF
0.13
0.13, 0.15, 0.12
0.23
0.23
Page on SITMUF
0.09
0.24, 0.63, 0.58
0.18
0.19
Combined Page MUF
0.10
0.52, 0.76, 0.79
0.21
0.21
Combined Page SITMUF
0.10
0.52, 0.76, 0.79
0.20
0.20
CUMUF
0.07
0.81, 0.05, 0.05
0.09
0.19
GEMUF
0.06
0.79, 0.41, 0.11
0.14
0.14
NP
0.15
0.98, 0.98, 0.98
0.44
0.62
SITMUF and CUMUF
0.06
0.79, 0.71, 0.75
0.18
0.13
SITMUF, Page, CUMUF
0.11
0.69, 0.57, 0.63
0.16
0.17
SITMUF, Page, GEMUF
0.11
0.52, 0.65, 0.79
0.19
0.21
Page, CUMUF, GEMUF
0.11
0.69, 0.57, 0.63
0.16
0.20Table 3
DPs forΣ3.
DP
Loss 1(0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3)
Loss 2(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3)
Loss 3(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)
Loss 4(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
MUF
0.07
0.66, 0.66, 0.66
0.15
0.17
SITMUF
0.05
0.81, 0.98, 0.81
0.43
0.20
Page on MUF
0.11
0.10, 0.10, 0.10
0.12
0.16
Page on SITMUF
0.07
0.20, 0.93, 0.20
0.83
0.18
Combined Page MUF
0.08
0.52, 0.98, 0.52
0.80
0.22
Combined Page SITMUF
0.07
0.52, 0.98, 0.52
0.79
0.22
CUMUF
0.06
0.80, 0.05, 0.80
0.05
0.07
GEMUF
0.06
0.91, 0.85, 0.96
0.18
0.20
NP
0.11
1.0, 1.0, 1.0
0.97
0.83
SITMUF and CUMUF
0.06
0.85, 0.97, 0.98
0.35
0.17
SITMUF, Page, CUMUF
0.10
0.24, 0.82, 0.88
0.14
0.16
SITMUF, Page, GEMUF
0.10
0.24, 0.82, 0.88
0.14
0.16
Page, CUMUF, GEMUF
0.10
0.09, 0.09, 0.09
0.11
0.15To illustrate the behavior of some of the tests usingΣ1 and loss 3 in Tables 1–3 ((0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)), Figure 5 plots the true loss and example MUF and SITMUF values (and the average SITMUF value over all simulations) in (a) and plots example DPs in (b). Figure 6 plots the residual standard deviation σ~j=σjj2-fΣ-1fT versus balance period for Σ1.Figure 5
ForΣ=Σ1, (a) the true loss, MB, SITMUF, and average SITMUF over all simulations versus balance period and (b) the alarm probability versus balance period (105 simulations) for SITMUF, Page1 SITMUF, CUMUF, and GEMUF.
(a)
(b)Figure 6
Residual standard deviationσ~i=σii2-fΣ-1fT versus balance period for Σ1.A sensitivity analysis was also performed by simulating 30% RSD inδ^R and 50% RSD in δ^S (see Section 3.2). With these relatively large RSDs in δ^R and δ^S, there is large uncertainty in the DP. For example, the 95% interval for DPs based on 105 simulations is (0.09, 0.41), (0.59, 0.99), (0.34, 1.0),(0.96, 1.0),(0.96, 1.0),(0.96, 1.0),(0.11, 0.92), (0.26, 0.97), (0.99, 1.0), (0, 56, 1.0), (0.35, 1.0), (0.93, 1.0), (0.31, 1.0), respectively, for the 13 tests for the loss in column 3 in Table 1 (0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0). The three least-sensitive DPs are in boldface and are for Page on SITMUF, combined Page on MUF, and combined Page on SITMUF. The most sensitive DP is also in boldface and underlined and is for CUMUF.


## 4.1. Sequential Statistical Testing of MB Sequences

Recall that NMA evaluates one or more MBs, where the MB is defined for balance periodj as MBj=(Ij-1+Tin,j-Tout,j)-Ij [38, 39]. Typically, many measurements are combined to estimate the terms Tin,Ibegin,Tout, and Iend in the MB; therefore, the central limit effect and years of experience suggest that MBs in most facilities will be approximately normally distributed with the mean equal to the true NM loss λj  and standard deviation σj, which is expressed as Xj~N(λj,σj), where X denotes the MB and the notation σj is a shortened version of σMB,j. A sequence of n MBs will be assumed to have approximately a multivariate normal distribution [38–43],  X1,X2,…,Xn~MVN(λ,Σ), where then-by-n covariance matrix is(12)Σ=σ12σ122⋯σ1n2σ212σ22⋯σ2n2⋮σn12σn22⋯σn2.The magnitude ofσj determines the amount of NM loss, λ=∑k=1nλk, which leads to high detection probability (DP). For example, suppose the facility tests for NM loss only, not for NM gain, and assume that Xj~N(λj,σj) is an adequate model. Then, if a false alarm probability (FAP) of α = 0.05 is desired, the alarm threshold is 1.65σj. In the case of testing for loss only, it follows that the loss detection probability 1-β for λ = 3.3σ and 1-β > 0.95 if λ> 3.3σj, where β is the nondetection (false negative) probability. The factor 3.3 arises from symmetry of the Gaussian distribution, requiring α = β = 0.05, and the fact that 1.65 = 3.3/2 is the 0.95 quantile of the N(0,1) distribution. One common goal is for DP = 1-β to be at least 0.95 if λ ≥ 1 SQ (significant quantity, which is 8 kg for Pu), which is accomplished if σj≤ SQ/3.3. If σj > SQ/3.3, this can be mitigated by reducing measurement errors to achieve σj≤ SQ/3.3 (if feasible) and/or by closing the balances more frequently, so there is less nuclear material transferred per balance period, which reduces σj [38, 39]. The DP of other safeguards measures such as enhanced containment and surveillance with smart cameras and/or remote radiation detection is difficult to quantify and is outside the scope of this paper. This section concludes with four remarks.Remark 1.
Large throughput facilities try to makeσj as small as reasonably possible and often try to keep σj small as a percent of throughput but cannot achieve σj≤ SQ/3.3. For example, suppose that there is a measurement error relative standard deviation of 0.5% of throughput. And suppose that the FAP/DP goals are α = 0.05 and 1-β = 0.95 and annual throughput is 100 SQ. Then, σj = 0.5 SQ = SQ/2 > SQ/3.3, so protracted diversion of 1 SQ over one year will not have a high DP. Therefore, one reason for frequent MB accounting is that abrupt diversion over hours or days is very likely to be detected [40]. As a complementary approach that is beyond the scope here, process monitoring [39] methods have the potential to detect off-normal facility operation that could misdirect NM to undeclared locations.Remark 2.
In the 1980s, some believed that a plant reporting a MB every 30 days would have higher DP than that same plant reporting a MB every year. However, [44] showed that, for optimal (from the diverter’s viewpoint) protracted diversion with the per-period loss being proportional to the row sums of the covariance matrix Σ of the MB sequence, annual MB testing has higher DP than monthly MB testing, and so, for such protracted diversion, “less frequent balance closure is better.” However, [44] conceded that NRTA has shorter detection times and higher DP against abrupt diversion. Publications [45–50] soon followed involving joint sequential tests: one tuned to detect protracted diversion and one more tuned for abrupt diversion. Such joint Page’s tests can be tuned to have high DP against abrupt loss while still having reasonably high DPs against protracted loss. Two types of combined Page’s tests are included in the simulation study in Section 4.3.Remark 3.
The assumptionX1,X2,…,Xn~MVN(λ,Σ) implies that Σ is known without estimation error. In practice, Σ is estimated using variance propagation applied to Xj=MBj=(Ij-1+Tin,j-Tout,j)-Ij and there will be estimation error in the estimate Σ^ [40].
The simulation study in Section4.3 includes a sensitivity study to assess the impact of estimation error Σ^ on the estimated false alarm probabilities (FAPs) and detection probabilities (DPs).Remark 4.
This section focuses on operator MB sequences. In international safeguards, the inspector randomly selects items to verify NM declarations made by the operator. The difference statistic,D=N∑j=1n((oj-ij)/n), defined as the average difference in the sample of sizen (extrapolated to the population of sizeN) between operator declarations (almost always based on operator measurements) and inspector measurements can be used as a test statistic, or the D statistic [2] can be used to compute the inspector MB statistic (or sequence). Inspector MB=MB-D, which could be analyzed using sequential statistical methods such as those in Section 4.3. Alternatively, one can test individually each of the n paired differences dj=oj-ij and the overall D statistic, and if none are found to be statistically significant, then the IAEA could rely on operator MB evaluation as in Section 4.3.


## 4.2. Propagation of Variance to EstimateΣ

EstimatingΣ is a key step required in frequent NMA. To illustrate, a simplified example model of a generic electrochemical facility with one input stream, one output stream, and one key inventory item will be used here [38]. Each measurement method is modeled here using a multiplicative measurement error model for the operator (O): Mi=μi(1+Si+Ri), with Si~N(0,δS2) and Ri~N(0,δR2), where Mi is the operator’s measured value of item i, μi is the true but unknown value of item i, Ri is a random error of item i, and Si is a short-term systematic error for item i. Then, the diagonal terms of Σ are calculated as(13)σi2=Tini2δin,R2+δin,S2+Touti2σout,R2+σout,S2+Ii2δinv,R2+Ii-12δinv,R2+Ii-Ii-12δinv,S2.The off-diagonal terms in Σ are calculated as(14)σij2=TiniTinjδin,S2+ToutiToutjδout,S2+IiIj+Ii-1Ij-1δinv,S2-IiIj-1δinv,S2+δinv,R2ifj-i=1-Ii-1Ijδinv,S2+δinv,R2ifi-j=1.In the last two terms, the random error of the inventory term is only applied if the condition is true. Reference [31] showed that the POV results for σi2 and σij2 are obtained by appropriate interpretation of GUM’s equation (1) in Section 2. For this simplified version of an example from [38], this leads to examples of 12-by-12 covariance matrices for monthly MBs over a one-year period; three example matrices Σ1, Σ2, and Σ3 are listed next.Three example covariance matrices,Σ1, Σ2, and Σ3, for a generic electrochemical facility, are given [38]:(15)1.00-0.480.010.010.01-0.481.00-0.480.010.010.01-0.481.00-0.480.010.010.01-0.481.00-0.480.010.010.01-0.481.00Equation (15) is Σ1, scaled to unit variance, only displaying 5 by 5 of the 12 by 12. This is the nominal case with 4 kg input per period and 40 kg inventory; δR = δS = 0.01.(16)1.000.170.330.330.330.171.000.170.330.330.330.171.000.170.330.330.330.171.000.170.330.330.330.171.00Equation (16) is Σ2, scaled to unit variance. This is the case with 4 kg input per period and 4 kg inventory.(17)1.000.580.670.670.670.581.000.580.670.670.670.581.000.580.670.670.670.581.000.580.670.670.670.581.00Equation (17) is Σ3, scaled to unit variance. This is the case with 4 kg input per period and 40 kg inventory and δS increased from 0.01 to 0.02.


## 4.3. Sequential Testing of Material Balance in Nuclear Material Accountancy

The assumptionX1,X2,…,Xn~MVN(λ,Σ) implies that Y=Σ-1/2X~MVN(Σ-1/2λ,I), where I is the identity matrix. The transform Y=Σ-1/2X is known in safeguards as the standardized independently transformed MUF (SITMUF, where MUF is another name for the MB), which is most conveniently computed using the Cholesky decomposition [43]. There are two main advantages of applying statistical tests to Y rather than to X. First, alarm thresholds depend only on the sequence length n forY and not on the form of the covariance matrix Σ. Because it is best to calculate thresholds using simulation, this is a logistic advantage. Second, the variance of the Y sequence decreases over time, so if a diversion occurs late in the analysis period, the DP is larger for the Y sequence than for the X sequence. Note that one cannot claim higher DP for the Y sequence than for the X sequence in general, because the true loss scenario is never known, and the DP can be larger for X than for Y for some loss scenarios, which is demonstrated in Section 4.The value ofYi can be calculated using Y=Σ-1/2X, but, more intuitively as the residual from the X sequence, Yj={Xj-E(Xj∣Xj-1,Xj-2,…,X1)}/σ~j, where the standard deviation σ~j is given by σ~j=σjj2-fΣ-1fT, where f=Σj,1:(j-1), the 1 to (j-1) entries in the jth row of Σ.Several reasonable statistical tests have been evaluated in [38, 41, 44–51] and are included in the simulation study in Section 4.4, including the following 13 tests:(1)
MUF test: this compares each MUF value to a threshold, which is the same as a Shewhart test in quality control (QC). The test alarms on periodj if MUFj/σj≥T for some threshold T(2)
SITMUF test: this compares each SITMUF valueY=Σ-1/2X to a threshold, which is the same as a Shewhart test in QC. The test alarms on period j if SITMUFj/σSITMUF,j≥T for some thresholdT(3)
Page’s test applied to MUF: Page’s test to test for loss is a sequence of sequential probability ratio tests, defined asPj=max⁡(0,Pj-1+xj/σj-k), where P0=0 [52]. The test alarms on period j if Pj>T. The parameter k is a control parameter that is optimal for detecting a shift from zero loss to loss λ if k=λ/2. The alarm thresholdT (usually denoted as h in literature on Page’s test, but this paper uses T for alarm threshold) is chosen so that the FAP per analysis period (usually one year) is 0.05 or whatever FAP is specified(4)
Page’s test applied to SITMUF: Page’s test to test for loss is a sequence of sequential probability ratio tests, asPj=max(0,Pj-1+yj-k). The alarm threshold T is chosen so that the FAP per analysis period (usually one year) is 0.05 or whatever FAP is specified(5)
Combined Page’s tests applied to MUF: the use of Page’s test with a large value ofk and small value of T has good DP for abrupt loss, and the use of Page’s test with a small value of k and large value of T has good DP for protracted loss. Therefore, a reasonable option is to use a combination of two Page’s tests, one with large k and one with smallk(6)
Applying combined Page’s tests to SITMUF(7)
CUMUF: at periodj, CUMUFj=∑i=1jxi is the sum of all MUF values from period 1 to j(8)
GEMUF: it has been shown that if the loss vectorλ is known, then the Neyman-Pearson test statistic is λTΣ-1x, which is known as a matched filter in some literature. The GEMUF statistic substitutes xT for λT, so GEMUF=xTΣ-1x. In simulation studies, λ is known, so the NP test statistic is useful for calculating the largest possible DP. The GE in GEMUF is a German abbreviation of Geschätzter, which means “estimated,” so GEMUF means estimated MUF, and GEMUF is the same as the Mahalanobis distance from the0 vector and Hotelling’s multivariate T statistic [51](9)
A nonsequential version of the Neyman-Pearson test,λTΣ-1x, is useful to calculate the largest possible DP for given Σ and λ. For completeness, four other combined tests are also considered(10)
SITMUF and CUMUF(11)
Page’s on SITMUF and CUMUF(12)
SITMUF, Page’s on SITMUF, and GEMUF(13)
Page’s on SITMUF, CUMUF, and GEMUFThis section concludes with four remarks.Remark 1 (SITMUF transform).
The SITMUF transform is recommended for two reasons. First, simulation is typically used to select alarm thresholds, and it is convenient to always work on the same scale when selecting alarm thresholds, so the fact thatY=Σ-1/2X~MVN(Σ-1/2λ,I) is convenient. Note that alarm thresholds could be selected on the basis of exact or approximate analytical results for some, but not all, of the tests. For example, there are approximate expressions for T andk [53]. Second, the standard deviation σ~j is given by σ~j=σj2-fΣ-1fT, where f=Σj,1:(j-1), the 1 to (j-1) entries in thejth row of Σ, so the standard deviation of the MUF residuals decreases in the later periods. Therefore, the independence transform is analogous to a bias adjustment, leading to smaller prediction variance in later periods, which tends to increase the DP for SITMUF compared to MUF (there are exceptions where the DP for MUF is larger than the DP for SITMUF; see Section 4.4, DP results).Remark 2 (choosing thresholds).
Thresholds can be chosen in many ways and can be assumed to be constant for each period or not. Therefore, simulation DP results in Section4.4 are not claimed to be optimal but are example DP results.Remark 3 (performance criteria).
The main performance criterion for comparing tests is the DP. But the average time to detection and robustness to misspecifying the covariance matrixΣ are also important.Remark 4.
There are other tests in the literature, including the power one and scan tests [38, 41–43].


## 4.4. Simulation Study

Example DP results are in Tables1–3 for Σ1, Σ2, and Σ3, respectively, using 105 simulations (so are repeatable to within ±0.01) to choose alarm thresholds and to estimate DPs. Tables 1–3 indicate that, as expected, there is no overall best test. However, Page’s test on SITMUF or a combined Page’s test on SITMUF is often among the best performers. Whether Page’s test on SITMUF has larger DP than Page’s test on MUF depends on the covariance matrix.Table 1
DPs forΣ1. Boldface entries have the largest DPs (excluding NP) for the respective column.
DP
Loss 1(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3)/10
Loss 2(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3)
Loss 3(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)
Loss 4(0,1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
MUF
0.07
0.58, 0.58, 0.58
0.15
0.19
SITMUF
0.20
0.65, 0.89, 0.81
0.82
0.58
Page on MUF
0.18
0.75, 0.85, 0.74
0.76
0.58
Page on SITMUF
0.71
0.82, 0.99, 0.56
1.0
0.99
Combined Page MUF
0.74
0.83, 0.99, 0.80
1.0
0.99
Combined Page, SITMUF
0.74
0.83, 0.98, 0.80
1.0
0.99
CUMUF
0.26
0.91, 0.64, 0.18
0.50
0.65
GEMUF
0.13
0.87, 0.70, 0.11
0.57
0.48
NP
0.82
0.99, 1.0, 0.99
1.0
1.0
SITMUF and CUMUF
0.23
0.80, 0.87, 0.80
0.81
0.62
SITMUF, Page, CUMUF
0.28
0.85, 0.85, 0.75
0.76
0.58
SITMUF, Page, GEMUF
0.20
0.84, 0.84, 0.72
0.74
0.61
Page, CUMUF, GEMUF
0.23
0.85, 0.75, 0.75
0.76
0.58Table 2
DPs forΣ2.
DP
Loss 1(0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3)
Loss 2(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3)
Loss 3(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)
Loss 4(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
MUF
0.07
0.60, 0.60, 0.60
0.18
0.18
SITMUF
0.06
0.62, 0.76, 0.80
0.14
0.14
Page on MUF
0.13
0.13, 0.15, 0.12
0.23
0.23
Page on SITMUF
0.09
0.24, 0.63, 0.58
0.18
0.19
Combined Page MUF
0.10
0.52, 0.76, 0.79
0.21
0.21
Combined Page SITMUF
0.10
0.52, 0.76, 0.79
0.20
0.20
CUMUF
0.07
0.81, 0.05, 0.05
0.09
0.19
GEMUF
0.06
0.79, 0.41, 0.11
0.14
0.14
NP
0.15
0.98, 0.98, 0.98
0.44
0.62
SITMUF and CUMUF
0.06
0.79, 0.71, 0.75
0.18
0.13
SITMUF, Page, CUMUF
0.11
0.69, 0.57, 0.63
0.16
0.17
SITMUF, Page, GEMUF
0.11
0.52, 0.65, 0.79
0.19
0.21
Page, CUMUF, GEMUF
0.11
0.69, 0.57, 0.63
0.16
0.20Table 3
DPs forΣ3.
DP
Loss 1(0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3)
Loss 2(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0)(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3)
Loss 3(0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)
Loss 4(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
MUF
0.07
0.66, 0.66, 0.66
0.15
0.17
SITMUF
0.05
0.81, 0.98, 0.81
0.43
0.20
Page on MUF
0.11
0.10, 0.10, 0.10
0.12
0.16
Page on SITMUF
0.07
0.20, 0.93, 0.20
0.83
0.18
Combined Page MUF
0.08
0.52, 0.98, 0.52
0.80
0.22
Combined Page SITMUF
0.07
0.52, 0.98, 0.52
0.79
0.22
CUMUF
0.06
0.80, 0.05, 0.80
0.05
0.07
GEMUF
0.06
0.91, 0.85, 0.96
0.18
0.20
NP
0.11
1.0, 1.0, 1.0
0.97
0.83
SITMUF and CUMUF
0.06
0.85, 0.97, 0.98
0.35
0.17
SITMUF, Page, CUMUF
0.10
0.24, 0.82, 0.88
0.14
0.16
SITMUF, Page, GEMUF
0.10
0.24, 0.82, 0.88
0.14
0.16
Page, CUMUF, GEMUF
0.10
0.09, 0.09, 0.09
0.11
0.15To illustrate the behavior of some of the tests usingΣ1 and loss 3 in Tables 1–3 ((0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0)), Figure 5 plots the true loss and example MUF and SITMUF values (and the average SITMUF value over all simulations) in (a) and plots example DPs in (b). Figure 6 plots the residual standard deviation σ~j=σjj2-fΣ-1fT versus balance period for Σ1.Figure 5
ForΣ=Σ1, (a) the true loss, MB, SITMUF, and average SITMUF over all simulations versus balance period and (b) the alarm probability versus balance period (105 simulations) for SITMUF, Page1 SITMUF, CUMUF, and GEMUF.
(a)
(b)Figure 6
Residual standard deviationσ~i=σii2-fΣ-1fT versus balance period for Σ1.A sensitivity analysis was also performed by simulating 30% RSD inδ^R and 50% RSD in δ^S (see Section 3.2). With these relatively large RSDs in δ^R and δ^S, there is large uncertainty in the DP. For example, the 95% interval for DPs based on 105 simulations is (0.09, 0.41), (0.59, 0.99), (0.34, 1.0),(0.96, 1.0),(0.96, 1.0),(0.96, 1.0),(0.11, 0.92), (0.26, 0.97), (0.99, 1.0), (0, 56, 1.0), (0.35, 1.0), (0.93, 1.0), (0.31, 1.0), respectively, for the 13 tests for the loss in column 3 in Table 1 (0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0). The three least-sensitive DPs are in boldface and are for Page on SITMUF, combined Page on MUF, and combined Page on SITMUF. The most sensitive DP is also in boldface and underlined and is for CUMUF.


## 5. Summary

Statistical analyses used to support safeguards conclusions require UQ, usually by estimating the RSD in random and systematic errors associated with each measurement method. This paper reviewed why UQ is needed in nuclear safeguards and examined recent efforts to improve both bottom-up and top-down UQ for calibration data. A key issue in bottom-up UQ using calibration with only a few calibration standards is that existing analytical approximations to estimate variance components are not sufficiently accurate, so this paper illustrated that simulation is needed. Once calibration UQ is well quantified, whenever improved bottom-up UQ predicts smaller measurement error RSDs than are observed in top-down UQ [1], this is evidence of significant unknown NDA error sources (“dark uncertainty,” [9]), which can then potentially be identified.The RSD of an assay method is often defined as the reproducibility standard deviation as estimated in an interlaboratory comparison. Recent options for top-down UQ, such as constrained estimators or Bayesian estimators that use prior information, offer possible improvements over existing variance component estimators. Any such improvements in estimated RSDs should be accompanied by uncertainties in the RSDs, which means that uncertainty in the estimated uncertainties matters [2, 7]. The ABC approach in Section 3 appears to provide a robust estimate of the posterior distribution of the RSDs.There are other types of bottom-up UQ used in safeguards not considered here. For example, the FRAM (fixed-energy, response function analysis with multiple efficiencies) gamma-based method [54] does not rely on calibration, and FRAM’s uncertainties are impacted by physical mismatch between test items and assay assumptions, which leads to item-specific bias, and also by uncertainties in nuclear data such as half-lives. This paper also used simulation to evaluate the impact of uncertainty in measurement error RSDs on estimated nuclear material loss detection probabilities in sequences of measured material balances. Many different sequential statistical tests were evaluated. In a related context, [2] evaluated the impact of uncertainty in measurement error RSDs on estimated DPs in verification data.



---

*Source: 2679243-2017-09-12.xml*