Monte Carlo Sampling Bias in the Microwave Uncertainty Framework

TL;DR

This paper investigates bias issues in the Microwave Uncertainty Framework's Monte Carlo sampling, proposing an alternative method that reduces bias at large sample sizes and analyzing its performance compared to the original.

Contribution

The paper identifies bias in the MUF's Monte Carlo sampling and introduces an alternative construction that asymptotically reduces bias for large sample sizes.

Findings

Bias is limited under certain conditions with the current method.

The proposed alternative reduces bias asymptotically at large sample sizes.

Neither method fully meets design goals at small sample sizes.

Abstract

Uncertainty propagation software can have unknown, inadvertent biases introduced by various means. This work is a case study in bias identification and reduction in one such software package, the Microwave Uncertainty Framework (MUF). The general purpose of the MUF is to provide automated multivariate statistical uncertainty propagation and analysis on a Monte Carlo (MC) basis. Combine is a key module in the MUF, responsible for merging data, raw or transformed, to accurately reflect the variability in the data and in its central tendency. In this work the performance of Combine's MC replicates is analytically compared against its stated design goals. An alternative construction is proposed for Combine's MC replicates and its performance is compared, too, against Combine's design goals. These comparisons are made within an archetypal two-stage scenario in which received data are first transformed in conjunction with shared systematic error and then combined to produce summary information. These comparisons reveal the limited conditions under which Combine's uncertainty results are unbiased and the extent of these biases when these conditions are not met. For small MC sample sizes neither construction, current or alternative, fully meets Combine's design goals, nor does either construction consistently outperform the other. However, for large MC sample sizes the bias in the proposed alternative construction is asymptotically zero, and this construction is recommended.

Equations

$${\vec{N}}_j^{\text{mbox}(\mathsf{T})}=F({\vec{Y}}_j,\vec{\nu })$$

$${\vec{M}}_{jq}^{\text{mbox}(\mathsf{T})}=F({\vec{Y}}_j,{\vec{S}}_q)$$

$$F(\vec{y},\vec{s})=\left(\begin{array}{c}f(y_1,s_1)\\ ⋮\\ f(y_K,s_K)\end{array}\right)$$

$${\vec{N}}^{\text{mbox}(\mathsf{C})}=\frac{1}{J}\sum _{j=1}^J{\vec{N}}_j^{\text{mbox}(\mathsf{T})}=\bar{F}({\vec{Y}}_{\bullet },\vec{\nu })$$

$${\vec{M}}_q^{\text{mbox}(\mathsf{C})}={\bar{M}}_{\bullet q}^{\text{mbox}(\mathsf{T})}+\frac{1}{\sqrt{J}}{\mathbf{U}}_{\text{mbox}\mathsf{C}}\sqrt{{\mathbf{D}}_{\text{mbox}\mathsf{C}}}{\vec{Z}}_q$$

$$\hat{\mbox{\boldmath$\Sigma$}}_{\vec{N}_{j}^{\mbox{\sf\tiny(T)}}}=\frac{1}{J-1}\sum_{j=1}^{J}(\vec{N}_{j}^{\mbox{\sf\tiny(T)}}-\bar{N}_{\bullet}^{\mbox{\sf\tiny(T)}})(\vec{N}_{j}^{\mbox{\sf\tiny(T)}}-\bar{N}_{\bullet}^{\mbox{\sf\tiny(T)}})^{\mathsf{T}}$$

$$\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}=\frac{1}{Q-1}\sum_{q=1}^{Q}(\vec{M}_{q}^{\mbox{\sf\tiny(C)}}-\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}})(\vec{M}_{q}^{\mbox{\sf\tiny(C)}}-\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}})^{\mathsf{T}}$$

$$E[{\bar{M}}_{\bullet }^{\text{mbox}(\mathsf{C})}]=E[\bar{F}({\vec{Y}}_{\bullet },{\vec{S}}_q)]$$

$$E[\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}]=V[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})].$$

$$E[\bar{F}({\vec{Y}}_{\bullet },{\vec{S}}_q)]=E[F({\vec{Y}}_j,{\vec{S}}_q)]$$

$$V[\bar{F}({\vec{Y}}_{\bullet },{\vec{S}}_q)]$$

$$E[{\bar{M}}_{\bullet }^{\text{mbox}(\mathsf{C})}]=E[\bar{F}({\vec{Y}}_{\bullet },{\vec{S}}_q)]$$

$$E[\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}]=V[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]+\frac{1}{J}\Psi$$

$$\Psi =V[F({\vec{Y}}_j,\vec{\nu })]-V[E[F({\vec{Y}}_j,{\vec{S}}_q)|{\vec{Y}}_j]].$$

$$E[{\bar{M}}_{\bullet }^{\text{mbox}(\mathsf{C})}]$$

$$\hat{\mbox{\boldmath$\Sigma$}}=\frac{1}{J-1}\sum_{j=1}^{J}(\vec{X}_{j}-\bar{X})(\vec{X}_{j}-\bar{X})^{\mathsf{T}}.$$

$$E[\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}]=V[\vec{M}_{q}^{\mbox{\sf\tiny(C)}}]-Cov[\vec{M}_{q}^{\mbox{\sf\tiny(C)}},\vec{M}_{q^{\prime}}^{\mbox{\sf\tiny(C)}}]$$

$$V[{\vec{M}}_q^{\text{mbox}(\mathsf{C})}]$$

$$E[\hat{}\mbox{\boldmath$\Sigma$}_{\vec{N}_{j}}^{\mbox{\sf\tiny(T)}}]=V[\vec{N}_{j}^{\mbox{\sf\tiny(T)}}]=V[F(\vec{Y}_{j},\vec{\nu})]$$

$$V[{\vec{M}}_q^{\text{mbox}(\mathsf{C})}]=V[\bar{F}({\vec{Y}}_{\bullet },{\vec{S}}_q)]+\frac{1}{J}V[F({\vec{Y}}_j,\vec{\nu })].$$

$$Cov[{\vec{M}}_q^{\text{mbox}(\mathsf{C})},{\vec{M}}_{q^{\prime }}^{\text{mbox}(\mathsf{C})}]$$

$$\mbox{relbias}[\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}]=\frac{E[\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}]-V[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]}{V[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]}$$

$$E[F({\vec{Y}}_j,{\vec{S}}_q)|{\vec{Y}}_j]$$

$$E[f(Y_{jk},S_{qk})|{\vec{Y}}_j]$$

$$E[\sin (Y_j+S_q)|Y_j=y]=\frac{\sin y\sin \delta }{\delta }.$$

$$V[E[\sin (Y_j+S_q)|Y_j]]=V\left[\frac{\sin Y_j\sin \delta }{\delta }\right]=\frac{{\sin }^2\delta }{{\delta }^2}V[\sin Y_j]$$

$$\Psi$$

$$\mbox{relbias}[\hat{\mbox{\boldmath$\Sigma$}}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}]=\frac{\frac{1}{J}\Psi}{V[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]}=\frac{\frac{1}{J}1}{\frac{1}{J}\frac{1}{2}+\frac{J-1}{J}0}=2.$$

$$k(y)=E[y^{S_q}]=\frac{y\sinh (\alpha \ln y)}{\alpha \ln y}.$$

$$Cov[F({\vec{Y}}_j,{\vec{S}}_q),F({\vec{Y}}_{j^{\prime }},{\vec{S}}_q)]=Cov[Y_j^{S_q},Y_{j^{\prime }}^{S_q}]$$

Full text

Keywords: Monte Carlo sampling, sampling bias, statistical software, statistical uncertainty propagation, systematic error, statistical experiment

1 Introduction

Modern national, industrial, and academic laboratories engaged in high-precision metrology rely on statistical software for multivariate and functional measurement uncertainty propagation and analysis. This software is typically highly complex and flexible and often has a Monte Carlo basis. Even in software well-designed from a statistical perspective, biases can be inadvertantly introduced due variously to flaws in statistical procedures, the algorithms that support them, or the algorithms’ coding. Statistical experiments are a natural and powerful way to test for such biases. We report the results of a case study of a microwave measurement uncertainty software package, called the Microwave Uncertainty Framework (MUF), in which a significant heretofore unknown bias in the software was detected, characterized, and corrected. This case study shows that elementary statistical performance testing can successfully identify such biases.

State-of-the-art microwave measurement relies on high-speed instrumentation including vector network analyzers (VNAs) operating in the frequency domain, temporal sampling oscilloscopes, and an array of other instruments, often used simultaneously in the same experiment. The refined measurements made possible by these arrangements allow investigators to, for example, identify the multiple reflections created by small imperfections in microwave systems, capture distortions due to the systems’ frequency-limited electronics, and study the role of noise. Statistical analysis of the data from this mix of instrumentation, including the conduct of uncertainty analyses, often involves shifts between the time and frequency domains. These shifts require that microwave uncertainty analyses account, particularly, for statistical correlations among the measurement uncertainties. To see this, consider that imperfections in microwave systems are often the source of unwanted reflections and attendant power losses. These temporal effects Fourier-map into the frequency domain as ripples with a characteristic period related to the inverse of the reflections’ time spacing. The VNA is currently the most accurate instrument for measuring these multiple reflections, and the errors made by this frequency sampling instrument typically manifest themselves as correlated time domain errors in the magnitudes, shapes, and positions of the multiple reflections. Statistical uncertainties in VNA measurements cannot be transformed correctly into the time domain without accounting for correlations created by the domain transformation [2].

Microwave measurement instrumentation with its often voluminous data production and the data-analytic need to track statistical correlations have motivated three automated approaches for statistical uncertainty analysis, the METAS VNA Tools II software package [3], the Garelli-Ferrero (GF) software package [4], and the MUF. The METAS and GF software packages support microwave multi-port network investigations, and offer fast, efficient sensitivity analysis implementations [5].

The MUF is a software suite created, supported, and made publicly available by the Radio Frequency Division of the U.S. National Institute of Standards and Technology. The MUF has capabilities similar to those of the METAS and GF software packages, while supporting radio frequency engineering applications beyond network analysis. The MUF’s general purpose is to provide automated multivariate statistical type A- and type B-evaluated [6] uncertainty propagation and analysis on a Monte Carlo (MC) basis [7] accessible through user-friendly interfaces. The MUF’s MC capability preserves non-Gaussian features of measured multivariate microwave signals, identifying systematic biases, for example, in signal calibration and processing steps. The MUF is composed of functional modules selectable by the user as needed for analyses. These include Model modules to flexibly represent microwave system elements. Model modules are useful, for example, for building calibration models, and they can be cascaded to represent increasingly complex systems. Other processing modules, termed here Transform modules, are available to perform oscilloscope and receiver calibrations, Fourier transforms, and other user-defined custom analytical transformations. Combine is a key module in the MUF, responsible for merging data, raw or transformed, to accurately reflect the variability in the data and in its central tendency. Combine is designed to be used at any point in extended analyses where repeated measurements must be merged. This flexibility is a powerful feature of the MUF.

This paper presents an analysis of the Combine module in the MUF. Because of the MUF’s distributed, multi-user, multi-purpose nature, Combine can be executed at different stages of uncertainty propagation analysis. We study a common use of Combine described by the two-stage scenario diagrammed in Fig. 1. In the first stage of this scenario, multivariate data are joined with shared systematic error in a bank of Transform modules and then, in the second stage, the transformed data are combined within the Combine module. Transform’s outputs take the form of nominal values of a selected mathematical transformation with associated uncertainties. Transform’s outputs in this form allow us in our two-stage scenario to study and assess Combine’s operation in which its inputs with their associated uncertainties are used to produce a summary mean output with an associated uncertainty. Combine represents the uncertainty in the summary mean in various fashions but provides the most detail in the form of a sample of MC replicates. Our analysis of Combine focuses specifically on the bias in the mean and covariance of these MC replicates. This analysis reveals that Combine’s construction of MC replicates is fundamentally biased, and we propose an alternative construction that effectively eliminates this bias.

The remainder of the paper is organized as follows. In Sect. 2 we analyze Combine’s performance in the two-stage scenario diagrammed in Fig. 1, showing that the sample mean of the MC replicates has zero bias and giving an analytical expression for the bias in the covariance of its MC replicates. This covariance bias is studied for specific cases of additive, multiplicative, exponential, and phase error. In Sect. 3 we propose an alternative construction for Combine’s MC replicates and show that the sample mean of Combine’s MC replicates has zero bias. We put tight bounds on the corresponding covariance bias and show that this bias is asymptotically zero; in this latter regard the proposed construction is better than the current method. Estimation bias is the primary concern in MC sampling, but MC estimation variability is also an issue. In Sect. 4 we continue our comparison of the current and alternative MC replicate constructions, comparing the variability in their sample means and sample covariances. We conclude in Sect. 5 with summary remarks supporting adoption of the proposed alternative MC replicate construction method in place of Combine’s current method. For the results presented in the following sections, we assume without note that the usual technical conditions pertain, that functions are measurable, that moments of sufficient order exist, etc.

2 Bias in the Combine module

We suppose in the two-stage scenario in Fig. 1 that the $J>1$ data vectors $\vec{Y}_{j}$ (of length $K$) are identically distributed and mutually independent and write \vec{Y}_{j}\sim(\vec{\mu},\mbox{\boldmath\Sigma}) to identify the mean $\vec{\mu}$ and covariance matrix $\Sigma$ of $\vec{Y}_{j}$. We also suppose that the MC-generated, length-$K$ errors $\vec{S}_{q}$, $q=1,...,Q$ are identically distributed, mutually independent, and independent of the sample of data vectors $\vec{Y}_{j}$. The mean and covariance of $\vec{S}_{q}$ are \vec{S}_{q}\sim(\vec{\nu},\mbox{\boldmath\Upsilon}). The covariance $\Sigma$ represents random uncertainty in the measurement of $\vec{\mu}$ while the errors $\vec{S}_{q}$ are systematic post-measurement errors introduced among the $\vec{Y}_{j}$ due, for example, to calibration adjustments.

Each data vector $\vec{Y}_{j}$ in Fig. 1 is operated on individually by Transform, producing for $\vec{Y}_{j}$ a nominal value

[TABLE]

for $j=1,...,J$ and a sample of vector Monte Carlo (MC) replicates

[TABLE]

for $q=1,...,Q$. The superscripts in (1) and (2) signify that these are Transform outputs in the first stage in our scenario. The MC replicates in (2) vary for a given $\vec{Y}_{j}$ only according to random replicates from the distribution of $\vec{S}_{q}$. The same $Q$ random replicates $\vec{S}_{q}$ are used to create each $\vec{Y}_{j}$’s sample of MC replicates. This models systematic errors that are shared among the $\vec{Y}_{j}$. The Transform module can similarly implement unshared systematic errors by using independent sets of $\vec{S}_{q}$ for each $\vec{Y}_{j}$, but the need for this capability rarely arises in application.

The transformation $F$ in Transform has the general form

[TABLE]

where $y_{k}$ and $s_{k}$ are the $k$th components of the vectors $\vec{y}$ and $\vec{s}$, respectively. The scalar-valued function $f$ is a user-specified parameter in Transform. Some choices of $f$ are $f(y,s)=y+s$, $f(y,s)=ys$, $f(y,s)=\sin(y+s)$, and $f(y,s)=y^{s}$, representing additive, multiplicative, phase, and exponential error, respectively. Transform also has many optional parameters, among them two matrix parameters, $\mathbf{T}_{Y}$ and $\mathbf{T}_{S}$. When, for example, the user specifies a matrix value for $\mathbf{T}_{Y}$, $F$ is applied to $\mathbf{T}_{Y}\vec{Y}$ instead of $\vec{Y}$. The matrix parameter $\mathbf{T}_{S}$ operates similarly. Mathematically, $\mathbf{T}_{Y}$ and $\mathbf{T}_{S}$ in the specification of $F$ are redundant; for example, we can without loss of generality take $\mathbf{T}_{Y}=\mathbf{I}$ by substituting $\mathbf{T}_{Y}\vec{\mu}$ and \mathbf{T}_{Y}\mbox{\boldmath\Sigma}\mathbf{T}_{Y}^{\mathsf{T}} for $\vec{\mu}$ and $\Sigma$, respectively. The optional use of $\mathbf{T}_{Y}$, though, gives Transform representational flexibility, allowing it, for example, to implement the Fourier transform of $\vec{Y}_{j}$. The covariance $\Sigma$ of $\vec{Y}_{j}$ is interpreted as the error covariance associated with measurement of $\vec{\mu}$, so $\Sigma$ and $\mathbf{T}_{Y}$ reflect unrelated physical processes and have distinct modeling roles. The matrices $\Upsilon$ and $\mathbf{T}_{S}$ play similarly distinct modeling roles.

Combine in Fig. 1’s two-stage scenario produces a nominal value $\vec{N}^{\mbox{\sf\tiny(C)}}$ for the transformed data and a sample of MC replicates $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ to describe the distribution and, particularly, the uncertainty of the transformed and combined data. These Combine outputs are given by

[TABLE]

and

[TABLE]

where $\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}=\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})$ and $\vec{Z}_{q}\sim\mbox{N}(\vec{0},\mathbf{I})$, $q=1,...,Q$. Further, the $\vec{Z}_{q}$ are independent of both the $\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}$ and $\mathbf{U}_{\mbox{\sf\tiny C}}\sqrt{\mathbf{D}_{\mbox{\sf\tiny C}}}$. The matrices $\mathbf{U}_{\mbox{\sf\tiny C}}$ and $\mathbf{D}_{\mbox{\sf\tiny C}}$ are the unitary and diagonal members, respectively, of the eigendecomposition \hat{\mbox{\boldmath\Sigma}}_{\vec{N}_{j}^{\mbox{\sf\tiny(T)}}}=\mathbf{U}_{\mbox{\sf\tiny C}}\mathbf{D}_{\mbox{\sf\tiny C}}\mathbf{U}_{\mbox{\sf\tiny C}}^{\mathsf{T}} of the sample covariance matrix

[TABLE]

associated with the nominal vectors $\vec{N}_{j}^{\mbox{\sf\tiny(T)}}$ at Combine’s input.

The vectors $\vec{Z}_{q}$ in (5) model standard normal variation along the principal axes of \hat{\mbox{\boldmath\Sigma}}_{\vec{N}_{j}^{\mbox{\sf\tiny(T)}}}. These standard normal variates are scaled by the standard deviations in $\sqrt{\mathbf{D}_{\mbox{\sf\tiny C}}}$ and then rotated by $\mathbf{U}_{\mbox{\sf\tiny C}}$ onto the coordinate axes of \hat{\mbox{\boldmath\Sigma}}_{\vec{N}_{j}^{\mbox{\sf\tiny(T)}}} to reflect the variability of the transformed data at Combine’s input. The components of the $\vec{Z}_{q}$ are chosen to be normally distributed based on the assumption that the number $J$ of Combine inputs and their independence are together great enough to support a Central Limit Theorem approximation. We present the $\vec{Z}_{q}$ as normally distributed because this is how they are generated in Combine. Only in subsection 4.2, however, is this distributional assumption necessary to our results.

The nominal value $\vec{N}^{\mbox{\sf\tiny(C)}}$ in (4) produced by Combine is a natural, intuitive summary of the central tendency of the transformed data provided that the transformed data are unimodal with little skew. The purpose of the MC replicates $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ is to indicate central tendency under more general conditions as well as to summarize the spread and distributional shape of the estimated central tendency. Formally, Combine is designed to produce a sample of MC replicates whose mean $\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}$ is an unbiased estimator of the vector $E[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]$ and whose covariance

[TABLE]

is an unbiased estimator of the covariance $V[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]$ of the vector $\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})$. In other words, the MC replicates in (5) should satisfy

[TABLE]

and

[TABLE]

We note for later use that under the conditions of our two-stage scenario the estimands in (8) and (9) can be expressed as

[TABLE]

and

[TABLE]

with $j\neq j^{\prime}$. Our analysis, summarized in Proposition 1 below, of the MC construction in (5) shows that Combine meets design goals (8) and (9) only under certain conditions, and that without these conditions Combine exhibits bias.

Proposition 1: Suppose that, in the two-stage scenario in Fig. 1 , we have $J>1$ independent, identically distributed data vectors \vec{Y}_{j}\sim(\vec{\mu},\mbox{\boldmath\Sigma}). Also suppose we have $Q$ independent, identically distributed errors \vec{S}_{q}\sim(\vec{\nu},\mbox{\boldmath\Upsilon}). Assume the sets of $\vec{Y}_{j}$ and $\vec{S}_{q}$ are independent. Suppose further that the Transform outputs $\vec{N}_{j}^{\mbox{\sf\tiny(T)}}=F(\vec{Y}_{j},\vec{\nu})$ and $\vec{M}_{jq}^{\mbox{\sf\tiny(T)}}$ are given by (1) and (2) with $F$ as in (3), and the Combine outputs $\vec{N}^{\mbox{\sf\tiny(C)}}$ and $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ are given by (4) and (5). Then

[TABLE]

and

[TABLE]

where $\Psi$ is the difference of two $K\!\times\!K$ covariances

[TABLE]

Proposition 1 establishes that the design goal in (8) is generally met by the MC replicates in (5), but the design goal in (9) is not. The covariance in the sample of MC replicates is biased by an amount $\Psi/J$. We will see in the next section that this bias can be positive or negative. We first prove the two parts (12) and (13) of the proposition.

Proof of (12): We first note that $E[\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}]=E[\vec{M}_{q}^{\mbox{\sf\tiny(C)}}]$ since the $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ are identically distributed. Then, conditioning on the factor $\mathbf{U}_{\mbox{\sf\tiny C}}\sqrt{\mathbf{D}_{\mbox{\sf\tiny C}}}$ in (5) and using that $\vec{Z}_{q}$ and $\mathbf{U}_{\mbox{\sf\tiny C}}\sqrt{\mathbf{D}_{\mbox{\sf\tiny C}}}$ are independent, we have

[TABLE]

Since $E[\vec{Z}_{q}]=\vec{0}$, this yields $E[\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}]=E[\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}]=E[\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})]$, which proves (12). ${\sqcap\sqcup}$

To prove (13) in the proposition, we need four lemmas, which we state here. Their proofs are given in the appendix. Lemma 1 concerns the sample covariance of cross-correlated vectors. Lemmas 2 and 3 are elementary conditioning argument-based results for auto- and cross-covariances. Lemma 4 is used here and in the proofs of subsequent propositions.

Lemma 1: Let \vec{X}_{j}\sim(\vec{\mu},\mbox{\boldmath\Sigma}) for $j=1,...,J$, $J>1$ with sample covariance

[TABLE]

Let $Cov[\vec{X}_{j},\vec{X}_{k}]=E[(\vec{X}_{j}-E[\vec{X}_{j}])(\vec{X}_{k}-E[\vec{X}_{k}])^{\mathsf{T}}]$ be the cross-covariance of $\vec{X}_{j}$ and $\vec{X}_{k}$, and suppose the $\vec{X}_{j}$ are cross-correlated with Cov[\vec{X}_{j},\vec{X}_{k}]=\mbox{\boldmath\Sigma}^{\prime} for all $j\neq k$. Then E[\hat{\mbox{\boldmath\Sigma}}]=\mbox{\boldmath\Sigma}-\mbox{\boldmath\Sigma}^{\prime}.

Lemma 2: Let $\vec{Z}\sim(\vec{0},\mathbf{I})$ be independent of the vector-matrix pair $(\vec{A},\mathbf{B})$. Then $V[\vec{A}+\mathbf{B}\vec{Z}]=V[\vec{A}]+E[\mathbf{B}\mathbf{B}^{\mathsf{T}}]$.

Lemma 3: Let $\vec{Z}_{1},\vec{Z}_{2}\sim(\vec{0},\mathbf{I})$, and suppose $\vec{Z}_{1}$, $\vec{Z}_{2}$, and $(\vec{A}_{1},\vec{A}_{2},\mathbf{B})$ are mutually independent. Then $Cov[\vec{A}_{1}+\mathbf{B}\vec{Z}_{1},\vec{A}_{2}+\mathbf{B}\vec{Z}_{2}]=Cov[\vec{A}_{1},\vec{A}_{2}]$.

Lemma 4: Let $\vec{S},\vec{S}^{\prime}$ be independent, identically distributed random vectors independent of the random vector $\vec{Y}$. Let $F(\vec{Y},\vec{S})$ be a vector function of $\vec{Y}$ and $\vec{S}$. Then $Cov[F(\vec{Y},\vec{S}),F(\vec{Y},\vec{S}^{\prime})]=V[E[F(\vec{Y},\vec{S})|\vec{Y}]]$.

Proof of (13): The MC replicate vectors $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ created by Combine are correlated with common cross-covariance $Cov[\vec{M}_{q}^{\mbox{\sf\tiny(C)}},\vec{M}_{q^{\prime}}^{\mbox{\sf\tiny(C)}}]$. Therefore, according to Lemma 1,

[TABLE]

with $q\neq q^{\prime}$. Using definition (5) for $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$, Lemma 2, definition (2) for $\vec{M}_{jq}^{\mbox{\sf\tiny(T)}}$, and the eigendecomposition \hat{\mbox{\boldmath\Sigma}}_{\vec{N}_{j}^{\mbox{\sf\tiny(T)}}}=\mathbf{U}_{\mbox{\sf\tiny C}}\mathbf{D}_{\mbox{\sf\tiny C}}\mathbf{U}_{\mbox{\sf\tiny C}}^{\mathsf{T}}, we have

[TABLE]

The Transform nominal values $\vec{N}_{j}^{\mbox{\sf\tiny(T)}}$ in (1) are independent and identically distributed so

[TABLE]

and (16) becomes

[TABLE]

Now consider the cross-covariance $Cov[\vec{M}_{q}^{\mbox{\sf\tiny(C)}},\vec{M}_{q^{\prime}}^{\mbox{\sf\tiny(C)}}]$ in (15). Using definition (5) for $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$, Lemma 3, and definition (2) for $\vec{M}_{jq}^{\mbox{\sf\tiny(T)}}$, we have

[TABLE]

the last equality holding because the data vectors $\vec{Y}_{j}$ are independent. Applying Lemma 4 to the covariance in (19) and substituting the result along with (18) back into (15) proves (13). ${\sqcap\sqcup}$

2.1 Example error models

Proposition 1’s point is that the MC replicates produced by Combine in our two-stage scenario have a covariance bias $\Psi/J$. In the remainder of this section we evaluate $\Psi$ for various error models, showing that $\Psi$ can be positive, negative, or zero. Where $\Psi$ is non-zero, we show in the univariate case $K=1$ that the relative bias

[TABLE]

approaches $\pm$20% in one example and even 200% in another.

Additive error: The function $f$ in (3) is $f(a,b)=a+b$ for additive error. In this case $F(\vec{y},\vec{s})=\vec{y}+\vec{s}$ and

[TABLE]

so $\Psi$ in (14) is identically zero. Thus for additive shared systematic error Combine’s MC replicates have both zero mean bias and zero covariance bias.

Multiplicative error: The function $f$ in (3) is $f(y,s)=ys$ for multiplicative error and the $k$th component of $F(\vec{Y}_{j},\vec{S}_{q})$ is $f({Y}_{jk},{S}_{qk})=Y_{jk}S_{qk}$ where $Y_{jk}$ and $S_{qk}$ are the $k$th components of $\vec{Y}_{j}$ and $\vec{S}_{q}$. We have

[TABLE]

Therefore $E[F(\vec{Y}_{j},\vec{S}_{q})|\vec{Y}_{j}]=F(\vec{Y}_{j},\vec{\nu})$ and $\Psi=0$. This shows that for multiplicative shared systematic error Combine’s MC replicates have both zero mean bias and zero covariance bias.

Phase error: The function $f$ in (3) is $f(y,s)=\sin(y+s)$ for phase error. In this case the covariance in Combine’s MC replicates can be biased. We focus on the univariate case $K=1$ in which we have scalars, $Y_{j}$ and $S_{q}$, and the phase error is uniformly distributed, $S_{q}\sim\mbox{Unif}(-\delta,\delta)$, $\delta>0$, with mean $\nu=0$ and range $2\delta$. We note first that

[TABLE]

Therefore

[TABLE]

and, using $\nu=E[S_{q}]=0$,

[TABLE]

To assess the relative size of the bias associated with $\Psi>0$ above, we consider the extremal case where $\delta=\pi$ and where $Y_{j}$ is $\pm\pi/2$ with equal probabilities. Then $\Psi=1$, $V[\sin(Y_{j}+S_{q})]=1/2$, and $Cov[\sin(Y_{j}+S_{q}),\sin(Y_{j^{\prime}}+S_{q})]=0$ for $j\neq j^{\prime}$. Using (LABEL:eq:simp2), we find that the relative bias (20) associated with the MC sample variance is

[TABLE]

Here the relative bias is 200% for any sample size $J$. This albeit exteme example demonstrates that very large relative biases are possible with Combine’s current method of MC replicate construction.

Exponential error: The function $f$ in (3) is $f(y,s)=y^{s}$ for exponential error. In this case the covariance in Combine’s MC replicates can be positively or negatively biased. We focus on the case $K=1$ of uniformly distributed scalars, $Y_{j}$ and $S_{q}$. For this case we find that $\Psi$ is broadly, but not always, negative.

Let $Y_{j}\sim\mbox{Unif}[a,b]$ and $S_{q}\sim\mbox{Unif}[1-\alpha,1+\alpha]$, $0\leq\alpha\leq 1$. We have $\nu=E[S_{q}]=1$ so $\Psi$ in (14) is $V[Y_{j}]-V[k(Y_{j})]$ where

[TABLE]

Using $V[Y_{j}]=\frac{1}{12}(b-a)^{2}$ and evaluating $V[k(Y_{j})]$ numerically, we find that $\Psi$ is slightly positive for small $b$, as shown in Fig. 2. Otherwise, in the region $0\leq a\leq b\leq 8$, $\Psi$ is negative, increasingly so for larger ranges $2\alpha$ and $b-a$.

In the cases presented in Fig. 2 for exponential error, the covariance

[TABLE]

is positive. According to (LABEL:eq:simp2), then, the relative bias (20) associated with $\Psi$ is strongest at the smallest sample size $J=2$, in which case

[TABLE]

Numerical evaluation of this expression yields the results presented in Fig. 3. At its strongest the relative bias approaches $\pm$20% for $\alpha=0.95$.

3 An alternative MC construction

The previous section shows that Combine’s MC replicates $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ in (5) generated for the two-stage scenario in Fig. 1 fail to fully meet Combine’s design goals (8) and (9). We propose in this section an alternative construction $\vec{M}_{q}^{\mbox{\sf\tiny(A)}}$ for Combine’s MC replicates. Like the $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ replicates in (5), the proposed $\vec{M}_{q}^{\mbox{\sf\tiny(A)}}$ replicates meet goal (8). Unlike the $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ replicates, the $\vec{M}_{q}^{\mbox{\sf\tiny(A)}}$ replicates essentially meet goal (9), doing so arbitrarily closely for sufficiently large MC replicate sample size $Q$.

Let

[TABLE]

where $\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}=\bar{F}(\vec{Y}_{\bullet},\vec{S}_{q})$ and $\vec{Z}_{q}\sim\mbox{N}(\vec{0},\mathbf{I})$, $q=1,...,Q$. Further, the $\vec{Z}_{q}$ are independent of both the $\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}$ and $\mathbf{U}_{\mbox{\sf\tiny A}}\sqrt{\mathbf{D}_{\mbox{\sf\tiny A}}}$. In this alternative construction the matrices $\mathbf{U}_{\mbox{\sf\tiny A}}$ and $\mathbf{D}_{\mbox{\sf\tiny A}}$ are now the unitary and diagonal members, respectively, of the eigendecomposition \hat{\mbox{\boldmath\Sigma}}_{\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}}=\mathbf{U}_{\mbox{\sf\tiny A}}\mathbf{D}_{\mbox{\sf\tiny A}}\mathbf{U}_{\mbox{\sf\tiny A}}^{\mathsf{T}} of the sample covariance

[TABLE]

associated with the means $\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}$ of the MC samples at Combine’s input. Proposition 2 below shows that basing the sample variability of the stage-two Combine MC replicates on the stage-one MC means $\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}$ instead of on the stage-one nominal values $\vec{N}_{j}^{\mbox{\sf\tiny(T)}}$ essentially removes the bias $\Psi/J$ identified in Proposition 1. This reduced bias is explained in some part by the greater information retained by using the MC means instead of the nominal values: the $M_{jq}^{\mbox{\sf\tiny(T)}}=F(\vec{Y}_{j},\vec{S}_{q})$ reflect nonlinearities in $F$ across the full distribution of $\vec{S}$, while the nominal values $\vec{N}_{j}^{\mbox{\sf\tiny(T)}}=F(\vec{Y}_{j},\vec{\nu})$ are only exposed to $F$ at the mean $\vec{\nu}$ of the $\vec{S}_{q}$ distribution.

Proposition 2: Let the set-up be the same as in Proposition 1 except that the Combine-stage MC replicates in Fig. 1 are given by $\vec{M}_{q}^{\mbox{\sf\tiny(A)}}$ in (21). Then

[TABLE]

and

[TABLE]

where $\Phi$ is the difference of two $K\!\times\!K$ covariances

[TABLE]

Proof: The proof of (23) is the same as that of (12) because $\vec{Z}_{q}$ and $\mathbf{U}_{\mbox{\sf\tiny A}}\sqrt{\mathbf{D}_{\mbox{\sf\tiny A}}}$ in (21) are again independent. To prove (24), we first note that the arguments based on Lemmas 1, 2, and 3 early in the proof of (13) apply also here, giving

[TABLE]

with

[TABLE]

and

[TABLE]

in which case

[TABLE]

Using Lemma 1, we write E[\hat{\mbox{\boldmath\Sigma}}_{\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}}] in (27) as

[TABLE]

Next,

[TABLE]

The $S_{q}$ in (31) are independent and identically distributed so, applying Lemma 4, $V[\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}]$ in (30) becomes, from (31),

[TABLE]

Similarly, the covariance $Cov[\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}},\bar{M}_{j^{\prime}\bullet}^{\mbox{\sf\tiny(T)}}]$ in (30) is

[TABLE]

Substituting (32) and (33) back into (30) yields

[TABLE]

Finally, substituting this back into (29) proves (24). ${\sqcap\sqcup}$

For the scalar case $K=1$ the relative bias

[TABLE]

associated with the MC sample variance in Proposition 2 has simple bounds, given in the following proposition. Following the proof of this proposition, we show by simple examples that these bounds are tight.

Proposition 3: The relative bias in (35) satisfies

[TABLE]

Proof: We prove first that the relative bias in (36) is non-negative. The variance of a random variable $X$ can be expressed by

[TABLE]

where $X,X^{\prime}$ are independent and identically distributed. Using conditional versions of (37), we have

[TABLE]

Define $\Delta(y,y^{\prime},s)=F(y,s)-F(y^{\prime},s)$. Then

[TABLE]

Variance is non-negative so $\Phi\geq 0$, proving the lower bound in Proposition 3. To prove the upper bound we note first that, applying Lemma 4 to the scalar case $V[\bar{F}(Y_{\bullet},S_{q})]$ of the target variance in (LABEL:eq:simp2), we have

[TABLE]

Then the relative bias given in (35) is

[TABLE]

establishing the upper bound in (36). ${\sqcap\sqcup}$

In the remainder of this section we look at the examples of additive and multiplicative error to see that the bounds in Proposition 3 on the relative bias of the MC sample variance \hat{\mbox{\boldmath\Sigma}}_{\vec{M}_{q}^{\mbox{\sf\tiny(A)}}} are tight.

Additive error: In this case $F(\vec{y},\vec{s})=\vec{y}+\vec{s}$ and $\Phi$ in (25) is

[TABLE]

This shows that Proposition 3’s lower bound is tight for the scalar case addressed there. More generally, it shows for additive shared systematic error that Combine MC replicates constructed according to (21) have both zero mean bias and zero covariance bias.

Multiplicative error: We focus on the univariate case $K=1$ in which we have scalars, $Y_{j}$ and $S_{q}$. Then

[TABLE]

The corresponding relative bias for our case $K=1$ is

[TABLE]

where in the last step we used the product rule for variance [9]. Expression (41) is less than or equal to $1/Q$, achieving $1/Q$ for $E[Y_{j}]=E[S_{q}]=0$, showing that the upper bound in Proposition 3 is tight. Thus, for multiplicative error with the alternative MC sample construction, the relative covariance bias can be as great as 100% in the extreme case $Q=1$; but for typical user choices of $Q$ it is no greater than a small fraction of a percent.

We conducted a computer experiment as a numerical check on the upper bound in Proposition 3 and, in particular, on (41) for relative bias in the case of multiplicative error. Combine MC samples with sizes ranging from $Q=3$ to $Q=300$ were constructed per the proposed alternative method for the univariate case ($K=1$) with independent standard normal datasets $Y_{j}$ of size $J=4$ and independent standard normal errors $S_{q}$. For each sample size 10,000 MC samples were created and their 10,000 sample variances were averaged. These averaged sample variances were compared against the target variance

[TABLE]

to estimate for each $Q$ the relative bias in the MC sample variance. The estimated relative biases are plotted in Fig. 4. According to (41), the estimated relative biases should agree with the solid red line in Fig. 4 given by $\mbox{relbias}=1/Q$. They do agree to within experimental uncertainty.

4 Relative variability

Sections 2 and 3 compared Combine’s current and proposed alternative constructions of MC replicates from the standpoint of bias in the MC sample means and covariances. Specifically, Propositions 1 and 2 established that the MC sample means $\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}$ and $\bar{M}_{\bullet}^{\mbox{\sf\tiny(A)}}$ are each unbiased. These propositions also established that MC sample covariance \hat{}\mbox{\boldmath\Sigma}_{\vec{M}_{q}^{\mbox{\sf\tiny(C)}}}^{2} is biased, while the bias of \hat{}\mbox{\boldmath\Sigma}_{\vec{M}_{q}^{\mbox{\sf\tiny(A)}}}^{2} is asymptotically zero as $Q\rightarrow\infty$. In this section we complete our comparison of the two constructions by considering the differences in the variabilities of their MC sample means and covariances.

4.1 Relative variability in MC sample means

The MC sample means with the current and alternative MC constructions (5) and (21) are, respectively,

[TABLE]

The difference in their degrees of variability (their covariances) is given by the following proposition.

Proposition 4: Consider the two-stage scenario in Fig. 1 with the same set-up as in Propositions 1 and 2. Then

[TABLE]

with the matrices $\Psi$ and $\Phi$ as defined in (14) and (25).

Proof: Lemma 2 yields for the MC sample means in (LABEL:eq:constructs) that

[TABLE]

in which case

[TABLE]

Using the definitions of $\Psi$ and $\Phi$ and the results in (17) and (LABEL:eq:resES2) for E[\hat{}\mbox{\boldmath\Sigma}_{\vec{N}_{j}^{\mbox{\sf\tiny(T)}}}] and E[\hat{}\mbox{\boldmath\Sigma}_{\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}}], we find

[TABLE]

Recalling definition (25) for $\Phi$, this proves the proposition. ${\sqcap\sqcup}$

For sufficiently large finite $Q$, the sign of the difference in the covariances $V[\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}]$ and $V[\bar{M}_{\bullet}^{\mbox{\sf\tiny(A)}}]$ is determined by $\Psi$ in (43). We saw in Sect. 2 that $\Psi$ can be positive, negative, or zero, so depending on the form of the Transform error model $F(\vec{Y}_{j},\vec{S}_{q})$ either of the two estimators $\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}$ or $\bar{M}_{\bullet}^{\mbox{\sf\tiny(A)}}$ can exhibit less variability. According to (43), when $\Psi$ is positive, the alternatively constructed MC replicates $\vec{M}_{q}^{\mbox{\sf\tiny(A)}}$ are better than the current $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ both because their sample mean $\bar{M}_{\bullet}^{\mbox{\sf\tiny(A)}}$ exhibits less variability and because their sample variance \hat{\mbox{\boldmath\Sigma}}_{\bar{M}_{j\bullet}^{\mbox{\sf\tiny(A)}}}^{2} is asymptotically unbiased. When $\Psi$ is negative, the picture is mixed: $\vec{M}_{q}^{\mbox{\sf\tiny(A)}}$ exhibits greater variability than $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$, but $\vec{M}_{q}^{\mbox{\sf\tiny(C)}}$ is biased. All these considerations of the relative variabilities of $\bar{M}_{\bullet}^{\mbox{\sf\tiny(C)}}$ and $\bar{M}_{\bullet}^{\mbox{\sf\tiny(A)}}$ are, of course, dominated by Proposition 4’s main import that the difference in their variabilties is asymptotically zero,

[TABLE]

4.2 Relative variability in MC sample covariances

We consider in this subsection the difference in the variabilities of the MC sample covariances with the current and proposed alternative constructions, limiting our considerations to the univariate ($K=1$) case. For $K=1$

[TABLE]

where $N_{j}^{\mbox{\sf\tiny(T)}}=F(Y_{j},\nu)$ and

[TABLE]

Even with the restriction $K=1$, assessing the variance of the sample variance of the $M_{q}^{\mbox{\sf\tiny(C)}}$ in (44) is difficult, necessarily involving fourth moments. Because the $Z_{q}$ in (44) are normal, the following lemma (proved in [8]) is useful.

Lemma 5: Let $X_{n}=\mu_{n}+\sigma Z_{n}$ for $n=1,...,N$ where $\sigma$ and the $\mu_{n}$ are constants and the $Z_{n}$ are mutually independent and standard normal-distributed. Let

[TABLE]

Then $E[S_{X}^{2}]=u^{2}+\sigma^{2}$ and

[TABLE]

Let $\bar{M}_{\bullet}^{\mbox{\sf\tiny(T)}}$ be the set of $Q$ MC sample means $\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}$. Conditioned on $\bar{M}_{\bullet}^{\mbox{\sf\tiny(T)}}$ and $S_{N_{j}^{\mbox{\sf\tiny(T)}}}^{2}$, Combine’s MC replicates $M_{q}^{\mbox{\sf\tiny(C)}}$ in (44) are independent and normally distributed with means $\bar{M}_{\bullet q}^{\mbox{\sf\tiny(T)}}$ and variance $\frac{1}{J}S_{N_{j}^{\mbox{\sf\tiny(T)}}}^{2}$. Then according to Lemma 5,

[TABLE]

and

[TABLE]

where

[TABLE]

Then

[TABLE]

A parallel calculation for the alternatively constructed MC replicates

[TABLE]

yields

[TABLE]

We have, therefore, from (46) and (47) that

[TABLE]

We now pursue expressions for the difference in the two variances in (48) in the cases of additive and multiplicative error.

Additive error model: For $f(y,s)=y+s$ we have $\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}=Y_{j}+\bar{S}_{\bullet}$ and $N_{j}^{\mbox{\sf\tiny(T)}}=Y_{j}+\nu$ where $\nu=E[S_{q}]$. Then $S_{\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}}^{2}=S_{Y_{j}}^{2}$ and $S_{N_{j}^{\mbox{\sf\tiny(T)}}}^{2}=S_{Y_{j}}^{2}$ so

[TABLE]

Multiplicative error model: For $f(y,s)=ys$ we have $\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}=Y_{j}\bar{S}_{\bullet}$ and $N_{j}^{\mbox{\sf\tiny(T)}}=Y_{j}\nu$. Then $S_{\bar{M}_{j\bullet}^{\mbox{\sf\tiny(T)}}}^{2}=\bar{S}_{\bullet}^{2}S_{Y_{j}}^{2}$ and $S_{N_{j}^{\mbox{\sf\tiny(T)}}}^{2}=\nu^{2}S_{Y_{j}}^{2}$. The factors $\bar{S}_{\bullet}^{2}$ and $S_{Y_{j}}^{2}$ are independent, so using the product rule for variance [9], we find

[TABLE]

Let $\tau$, $\omega$, and $\psi$ be the second, third, and fourth central moments of $S_{q}$, and let $\phi$ be the fourth central moment of $Y_{j}$. Moment results for the sample mean and sample variance in [10, 11] then yield

[TABLE]

It then follows from (50) that

[TABLE]

This suggests generally for multiplicative error that, for $Q\rightarrow\infty$,

[TABLE]

Computer experiments were done to check the approximation in (48) for the difference in the variances of the MC sample variances. According to results (49) and (52), the difference in the variances should be zero and asymptotically zero for additive and multiplicative noise, respectively. These results were confirmed in experiments with different distributions for $Y_{j}$ and $S_{q}$ and different data sample sizes $J$. The results in Fig. 5 are typical. The plots in Fig. 5 show the estimated relative difference

[TABLE]

the top plot for a case of additive noise and the bottom plot for a case of multiplicative noise. In each plot the data and noise are standard normal-distributed and $J=4$. The plotted relative differences defined by (53) have a potential range of $\pm 100\%$. Each plotted point in Fig. 5 is estimated from an independent set of 100,000 computer trials. The experiment results in Fig. 5 and corresponding results obtained for other distributions and sample sizes confirm (49) and (52) under broad conditions.

Our aim with results (49) and (52) was to discover whether either method of constucting MC replicates dominates the other with respect to variance of the MC sample variance. These results, and the experiments confirming them, indicate that neither MC replicate construction method dominates the other when the error is additive or multiplicative. For error models beyond the additive and multiplicative cases, little analytical headway seems possible, so we turn directly to computer experiments.

Presented in Fig. 6 are experiment results obtained for the phase and exponential error models. The top plot in Fig. 6 shows the relative difference (53) of variances for phase error for the extremal case discussed in Sect. 2 in which the error $S_{q}$ is distributed $\mbox{Unif}[-\pi,\pi]$, the data $Y_{j}$ are $-\pi/2$ and $\pi/2$ with equal probabilities, and $J=4$. The middle and bottom plots show two cases of exponential error, both in which the data and exponential error are uniformly distributed, $Y_{j}\sim\mbox{Unif}[0,b]$ and $S_{q}\sim\mbox{Unif}[1-\alpha,1+\alpha]$, with $b=8$ and $\alpha=.95$ in the middle plot and $b=1$ and $\alpha=.95$ in the bottom plot. These two cases of exponential error are the two cases in Fig.  3 where the relative bias of Combine’s MC sample variance is most extreme—approaching $\pm$20%.

The results in Fig. 6 show that neither MC construction method dominates the other by having consistently smaller variance in its sample variance. In the middle plot the relative difference (53) in variances is negative, meaning that the sample variance with the current method has less variability. In the bottom plot, though, also with the exponential error model, the sample variance with the alternative method has less variability. Also, the three plots illustrate that the relative difference in variances can exhibit different degrees and types of transient behavior for small $Q$. The bottom plot shows almost no transient change, while the middle plot shows significant transient change before settling toward a limiting non-zero relative difference. The top plot shows that the relative difference in the variances can even change sign before approaching its limit value.

5 Summary remarks

The MUF is a powerful tool for uncertainty modeling and analysis relating to data obtained in high-precision microwave experiments, and the Combine module is a key component of the MUF. We compared the MC replicates currently constructed by Combine with those based on an alternative construction, using bias and variance of the MC sample mean and sample covariance as performance measures. We showed first that, with the current method of Combine MC replicate construction, the MC sample covariance is biased. Examples showed that this bias can be unacceptably large—200% in one extreme example and approaching $\pm$20% in others—and cannot be reduced to a tolerable level by choosing the MC sample size $Q$ sufficiently large. The MC sample covariance using the alternative construction for MC replicates is also biased, but this bias is asymptotically zero with $Q$.

Bias is the primary concern in MC sampling and the distinction, the current method being biased and the alternative being asymptotically unbiased, is the two construction methods’ most important difference. Looking beyond bias to the difference in the variabilities of the MC sample means with the two methods, we showed that this difference is asymptotically zero, with neither method dominating the other for small $Q$. Comparing the variabilities of the MC sample variances was similarly nuanced and non-determinative: in the cases of additive and multiplicative error, the difference in the variances of the sample variances is zero or asymptotically zero, while for phase and exponential error, neither method consistently out-performs the other.

We showed in this case study of bias in uncertainty propogation software that our proposed alternative MC replicate construction method has an important advantage with regard to MC sample covariance bias, while lacking any clear disadvantage relative to the current method. Consequently, the current method of constructing MC replicates in the MUF Combine module is set to be replaced with our proposed alternative. This study shows both that unknown, inadvertent biases are potentially present in even well-designed statistical software and that statistical experiments can successfully identify these biases. We urge that statistical performance tests be standard for modern software uncertainty propagation tools, and we anticipate that the statistical approach used here will be useful to future analyses of MUF performance and of the performance of other similar statistical software for uncertainty propagation.

Acknowledgment

This investigation benefited from early-stage discussions with Chih-Ming Wang and Sarah Streett, members of the Statistical Engineering Division of the National Institute of Standards and Technology.

Appendix

We prove the four lemmas used in the proof of Proposition 1.

Lemma 1: Let \vec{X}_{j}\sim(\vec{\mu},\mbox{\boldmath\Sigma}) for $j=1,...,J$, $J>1$ with common cross-covariance Cov[\vec{X}_{j},\vec{X}_{k}]=\mbox{\boldmath\Sigma}^{\prime} for all $j\neq k$. Then E[\hat{\mbox{\boldmath\Sigma}}]=\mbox{\boldmath\Sigma}-\mbox{\boldmath\Sigma}^{\prime} where \hat{\mbox{\boldmath\Sigma}} is the sample covariance matrix

[TABLE]

Proof:

[TABLE]

Lemma 2: Let $\vec{Z}\sim(\vec{0},\mathbf{I})$ be independent of the vector-matrix pair $(\vec{A},\mathbf{B})$. Then $V[\vec{A}+\mathbf{B}\vec{Z}]=V[\vec{A}]+E[\mathbf{B}\mathbf{B}^{\mathsf{T}}]$.

Proof:

[TABLE]

Lemma 3: Suppose $\vec{Z}_{1}$, $\vec{Z}_{2}$, and $(\vec{A}_{1},\vec{A}_{2},\mathbf{B})$ are mutually independent with $\vec{Z}_{1},\vec{Z}_{2}\sim(\vec{0},\mathbf{I})$. Then $Cov[\vec{A}_{1}+\mathbf{B}\vec{Z}_{1},\vec{A}_{2}+\mathbf{B}\vec{Z}_{2}]=Cov[\vec{A}_{1},\vec{A}_{2}]$.

Proof:

[TABLE]

Lemma 4: Let $\vec{S},\vec{S}^{\prime}$ be independent, identically distributed random vectors independent of the random vector $\vec{Y}$. Let $F(\vec{Y},\vec{S})$ be a vector function of $\vec{Y}$ and $\vec{S}$. Then $Cov[F(\vec{Y},\vec{S}),F(\vec{Y},\vec{S}^{\prime})]=V[E[F(\vec{Y},\vec{S})|\vec{Y}]]$.

Proof:

[TABLE]

The two functions $E[F(\vec{Y},\vec{S})|\vec{Y}=\vec{y}]$ and $E[F(\vec{Y},\vec{S}^{\prime})|\vec{Y}=\vec{y}]$ are identical so the first covariance on the right in (5) is $V[E[F(\vec{Y},\vec{S})|\vec{Y}]]$. Also, $F(\vec{Y},\vec{S})$ and $F(\vec{Y},\vec{S}^{\prime})$ are conditionally independent given $\vec{Y}$, so their conditional covariance on the right in (5) is zero. ${\sqcap\sqcup}$