Efficiency for , for selected polynomial orders and numbers of parameters.
Statistical signal processing  traditionally focuses on extraction of information from noisy measurements. Typically, parameters or states are estimated by various filtering operations. Here, the quality of signal processing operations will be assessed by evaluating the statistical uncertainty of the result . The processing could for instance simulate, correct, modulate, evaluate, or control the response of a physical system. Depending on the addressed task and the system, this can often be formulated in terms of a differential or difference signal processing model equation in time, with uncertain parameters and driven by an exciting input signal corrupted by noise. The quantity of primary interest may not be the output signal but can be extracted from it. If this uncertain dynamic model is linear-in-response it can be translated into a linear digital filter for highly efficient and standardized evaluation . A statistical model of the parameters describing to which degree the dynamic model is known and accurate will be assumed given, instead of being the target of investigation as in system identification . Model uncertainty (of parameters) is then propagated to model-ing uncertainty (of the result). The two are to be clearly distinguished – the former relate to the input while the latter relate to the output of the model.
Quantification of uncertainty of complex computations is an emerging topic, driven by the general need for quality assessment and rapid development of modern computers. Applications include e.g. various mechanical and electrical applications [5-7] using uncertain differential equations, and statistical signal processing. The so-called brute force Monte Carlo method [8-9] is the indisputable reference method to propagate model uncertainty. Its main disadvantage is its slow convergence, or requirement of using many samples of the model (large ensembles). Thus, it cannot be used for demanding complex models. The ensemble size is a key aspect which motivates deterministic sampling. Small ensembles are found by substituting the random generator with a customized deterministic sampling rule. Since any computerized random generator produces a pseudo-random rather than a truly random sequence, this is equivalent of modifying the random generator to be accurate for small ensembles of definite size, rather than being asymptotically exact (infinite ensembles). Correctness of very large ensembles is of theoretical but hardly practical interest for complex models, if the convergence to the asymptotic result is very slow.
2. Modeling uncertainty of signals
2.1. Problem definition
Suppose the (output) signal of interest is generated from the (input) signal passing through a dynamic system , with parameters ,
The model is given in uncertain parameters, which can be arranged in a column vector . For systems continuous-in-time (CT), is the differential operator in time while for systems discrete-in-time (DT), is the negative unit displacement operator, . There are several approximate methods to sample CT systems to DT systems, see  and references therein. The discretization techniques are beyond the scope of this presentation and DT systems will be assumed. If , there is feedback in the system which results in an impulse response of infinite duration. For finite accuracy however, the duration is finite. The system is linear-in-response, . Most importantly, the system is non-linear-in-parameters if . This is the typical situation addressed here.
Systems of the form in Eq. 1 may be directly realized as digital filters,, where denotes the filtering operation. The coefficients and are the numerator and denominator coefficients of the filter with impulse response , respectively. Its z-transform is obtained with the substitution . The parameterization can be changed to for instance gain, poles and zeros , or poles and residues ,
The parameterization should be carefully chosen as it affects the convergence rate of Taylor expansions (section 3.1) as well as the physical interpretation. The parameters and their statistics are preferably extracted from measurements using system identification techniques . Note that complex-valued poles and zeros are conjugated in pairs .
The problem to be addressed is the statistical evaluation of any function , given statistical models of and . It will here consist of evaluating its time-dependent mean and standard deviation . Without loss of generality, the analysis will be made for . Digital filtering will be utilized for evaluating samples of the model, i.e. filtering with definite sets of and signals .
Statistical expectations of any signal, model or function over finite discrete as well as continuous ensembles or probability distributions (no subscript) are defined as,
Samples of are labeled , with their components organized in columns. Sample indices will be given as superscripts in parenthesis, eg. is a column vector denoting the sample of parameter . Variations from the mean are written as .
Only uniform and normal distributions will be utilized. Either the mean and standard deviation, or the interval in brackets will be given in parenthesis, e.g. . Statistical moments carry the information contained in the marginalized probability density functions (pdf) , where denotes the sample space. While describes the width of , is related but different to its skewness . Further, the shape is reflected in , similarly to the curtosis . Since and are normalized and symmetric , and . Their differences are first reflected in their fourth moment, for and , respectively. The maximum variation of the parameter is expressed by the range . Dependencies are expressed in mixed moments . The discussion will be limited to correlations described by the covariance matrix , where the vector multiplication is an outer product.
Matrix size will be indicated with subscripts, e.g. is a matrix of rows and columns with elements , and . The identity matrix will be denoted, while matrices with equal elements will have their size attached, . For a matrix (vector) , is a vector (diagonal matrix) with components (diagonal elements) equal to the diagonal elements (components) of . The trace of a matrix is denoted .
A method will be stated intrusive if manipulations of the model are required. For the targeted highly complex models, it will be assumed that the computational cost for their evaluation dominates all other calculations. The efficiency of any method will accordingly be defined by the least required number of evaluations of the original model.
2.3. Fundamentals of non-linear propagation of uncertainty
Linearity in parameters (LP) is to be distinguished from linearity in response (LR),
for some vector . Different concepts of linearity are used, for LP models. Strictly speaking, LP denotes models that are affine, i.e. written as linear combinations of their parameters. Most constructed systems are designed to be as close to LR as possible while most models are not LP. There is hence no contradiction in non-linear (LP) propagation of uncertainty with linear (LR) digital filters, as here.
For non-linear propagation of uncertainty, the asymmetry of the resulting pdf is central. It can be expressed as a lack of commutation of non-linear propagation and statistical evaluation of a center value , as measured with the scent ,
The method for evaluating the center is left unspecified, as there are several alternatives. The most common choice is to use the mean, . The lowest order approximation of the scent can then be obtained by calculating the expectation of a Taylor expansion (section 3.1), , where is the Hessian matrix signal of , evaluated at . The scent is related to the skewness . The additional asymmetry caused by the non-linearity of the model is measured with the scent but differently. The scent addresses how parametric uncertainties are propagated and not how the result is distributed, e.g. for all LP models for which may attain any value. A finite scent thus implies the model is not LP, but not the reverse. The scent should not be confused with bias. Bias is a property of an estimator, while scent is a property of a model. For every model, such as the REF (section 6.1), many different estimators of can be used, e.g. the different ensembles in section 5.6, see result in Fig. 5 (left). Consequently, an unbiased estimator of correctly accounts for rather than ignores its finite scent, or deviation from .
The scent is important since and not is the main result utilized in applications. The corresponding difference  in the standard deviation from its linearized approximation , with , affects the confidence in the result. Its accuracy is usually less critical. An accurate evaluation of the scent is perhaps the strongest feature of the unscented Kalman filter, which provides the foundation for the presented approach as well as the origin of the term ‘scent’.
3. Conventional methods
A brief resume of the most traditional related methods of uncertainty propagation, applicable to signal processing models, is here given together with their pros and cons. Advanced intrusive methods like e.g. polynomial chaos expansions [14-15] not directly related to the proposed method are omitted.
3.1. Taylor expansions
The indisputable default methods of uncertainty propagation are based on Taylor expansions. These methods are intrusive if the differentiations are made analytically. Convergent series require regular differentiable models and numerical or analytical complexity make them error prone. Their applicability is therefore limited for complex models.
The transfer function of the digital filter can be expanded in a Taylor series,
This defines sensitivity systems (column vector) , unique quadratic variation systems (matrix) , and so on. These variation systems differ (intrusive) from but may nevertheless be realized as digital filters [3,7,10], just as . The corresponding variation of is given by,
where are the impulse responses of the systems . Utilizing digital filters with impulse responses , the differentiations are conveniently done once, and not repeatedly for every signal . The linearity in parameters of the model can easily be studied for many different input signals , by evaluating . Due to the large number of variation systems, higher order perturbation analyses rapidly become intractable though. The established method is limited to linearization (LIN)  . It will always incorrectly yield vanishing scent, . A first order estimate of is instead given by the expectation of second term in Eq. 7, , where the matrix of Hessian signals is obtained with repeated digital filtering.
3.2. Brute force Monte Carlo
Monte Carlo (MC) methods [8-9], or random sampling of uncertain models was originally introduced and phrased ‘statistical sampling’ by Enrico Fermi already in the 1930’s . The MC methods realize uncertain signal processing models in finite ensembles. Every ensemble consists of a possible set of well-defined model systems, all (usually) having the same structure but slightly different parameter values. In the original so-called brute force Monte Carlo method, each set of parameters is assigned to the output of random generators with appropriate statistics. The convergence to the assigned statistics is very slow  but it is asymptotically exact and the required number of samples is essentially independent of the number of parameters. Hence it does not suffer from the curse-of-dimensionality of many other methods. The outstanding simplicity in application is likely the cause of its popularity, just as the slow convergence or low efficiency is the main reason for its failures.
In MC, arbitrary distributions and dependencies are usually obtained by means of transformations of samples of elementary distributions. Independent samples of any probability density function (pdf) can be constructed with the inverse transform method . It consists of a calculation of the inverse of its cumulative distribution function (cdf) and generation of a uniformly distributed random sequence ,
Covariance may be included with an appropriate transformation of samples of canonical parameters with ,
The matrices are found by calculating the eigenvalues and eigenvectors  of . This transformation makes the marginal pdfs to differ substantially from the univariate pdfs of the independent but scaled parameters ,
All are hence mixed according to . Dependencies are thus difficult to account for. One rare exception is provided by the multinomial distribution . It is often better to assign the pdfs to the canonical parameters in the original instead of the canonical basis. The transformation then reads . As required, it leaves invariant. The marginalization in Eq. 10 changes accordingly, . Since the transformation of contains cancelling operations and , it is generally less distorting than . Indeed, if the commutator vanishes, . The transformation must satisfy the stronger criterion to avoid mixing. For any transformation , an indicator of mixing of the components of is given by,
A simple example illustrates that the mixing effect can be considerable, even for minute correlations. Assume a model has two parameters with a covariance matrix,
Large rotations are required because the canonical variances are similar, i.e. is almost degenerate. As shown in Fig. 1, the large rotations mix the assigned pdfs to marginal pdfs beyond recognition for the transformation but not for .
Specifying both marginal probability distributions and covariance is either redundant or inconsistent, as the latter is uniquely determined by the former. Nevertheless, this reflects the typical available information for signal processing applications. The moments can be accurately determined  for sufficiently large data sets but the joint distribution is hardly ever known with any precision. Some of its properties are usually assigned, with varying degree of confidence. For instance, the allowed maximal range of the parameters of digital filters is given by stability constraints. The transformation technique above is well adapted to these facts, since the covariance is prioritised. The transformation will be utilized in section 5.2 to include correlations with limited mixing of the statistics assigned to independent normalized canonical parameters .
3.3. Refinements of Monte Carlo
To increase the efficiency of MC, the original brute force sampling technique has been further developed in mainly two directions: model simplification and sample distribution improvement. In response surface methodology (RSM) , the model is replaced by a simple approximate surrogate model. A model of order may be found by applying linear (with respect to ) regression at collocation points  ,
where is basis function . Since it may be non-linear, RSM allows for non-linear propagation of uncertainty and may give a substantially different and more accurate result than LIN. If only linear basis functions are used , RSM becomes equivalent to LIN. The best least square approximation is directly obtained from Eq. 13 ,
Let utilize a complete set of mixed polynomial basis functions up to order . Its least number of collocation points grows rapidly with both the number of parameters and polynomial order ,
In practice, often yields an unacceptable number of samples, see table 1.
The distribution of samples may be improved with stratification, as in Latin Hypercube sampling (LHS) . By dividing the sample space into intervals, or stratas representing equal probability the need for large ensembles is reduced. In LHS, each parameter is sampled exactly once in each of its stratas giving a generalized latin square . This selection pushes the samples away from each other and distributes them more evenly. To illustrate the improvement with stratification, sample one parameter . After division into intervals of equal probability, samples are found with the inverse transform method described in section 3.2 (Eq. 8). As seen in Fig. 2, the convergence improves dramatically. Still, even for samples the second moment (left) varies noticeably. The convergence is generally poorer for higher order moments , as shown for (right).
In this case, it is questionable if 100 samples are sufficient to represent as few as four moments . The probabilistically evenly distributed fixed grid (dashed) converges more rapidly to the proper statistics. Despite the prevailing tradition, there is no absolute requirement of using a random generator to represent statistical information. Fixed grids are examples of deterministic sampling. Stratification provides an interesting intermediate type of sampling since it is partially deterministic – the strata are constructed deterministically but the samples within each stratum are generated randomly. The construction of a fixed grid requires focus on the most relevant features. To reproduce exactly, a very sparse grid or few deterministic samples are needed,
If the problem at hand only depends on these moments, the exact solution will be obtained. The size of such small ensembles must be fixed, no matter how they are generated. Adding, or perturbing a single sample would modify the statistics substantially.
4. Deterministic sampling
Deterministic sampling (DS) of uncertain systems is a viable alternative to random sampling (RS). Instead of using random generators, specific DS rules are devised to generate appropriate, but still statistical (Fermi’s notation, see section 3.2) ensembles. A rudimentary example illustrates the principle: Assume a model depends on one parameter with mean and variance . To estimate the mean and the variance of the model, the samples (filter parameters) are appropriate since they satisfy the desired statistics, and . The formula for constitutes the sampling rule and is the statistical ensemble containing only two model samples. By paying the computational cost of using more samples and improving the sampling rule, additional moments or other statistical features can be accounted for.
In deterministic sampling the model evaluations involve no approximations and are non-invasive. In many respects, deterministic sampling is constructed and optimized for quantification of modeling uncertainty: Minimal ensembles allow for evaluation of the most numerically demanding models. The model evaluations are exact and non-invasive to fully respect non-linear deeply hidden parameter dependences. Only vaguely known statistics of the model is approximated.
4.1. Concepts of deterministic sampling
DS does not per se specify the goal of sampling, e.g. given mean and covariance of the parameters. In the example at the end of section 3.3, the primary target was the joint pdf of the parameters. In section 4.2, the target is . In section 5, this will be complemented with additional requirements. DS can also be utilized for direct evaluation of confidence intervals . The targets of various DS methods may differ but the focus on the most influential statistical aspect and customization is shared. In stark contrast, almost without exception RS targets the joint pdf of the parameters and ignores the final utilization. Adaptation and fixed ensemble sizes provides the principal means to improve the efficiency of sampling.
4.2. Propagation of covariance in the standard unscented Kalman filter
The reference will be the specific variant of DS used for propagating covariance in what will be referred to as the standard unscented Kalman filter (UKF) [21-23]. The ensemble consists of samples, or sigma-points,
where denotes the -th column of . The sampling rule is manifested in the square root calculation of the covariance matrix . As suggested  it may be found with a Cholesky factorization . The square root matrix is not unique though – the Cholesky root is upper triangular and thus asymmetric. A more symmetric standard alternative is to evaluate the matrix square root in a canonical basis  where is diagonal. The canonical variations will be unit vectors in the positive and negative directions of the principal axes of the covariance matrix, amplified by the marginal standard deviations and most importantly, . For many parameters with large covariance, the scaling with may cause the UKF to fail since the scaling is not related to the variability of the parameters, only their total number. A possible solution to the scaling problem is provided by the scaled unscented transformation . However, it is based on Taylor expansions and thus suffers from an approximation problem of the model.
5. Sampling with conservation of moments
One class of methods of deterministic sampling conserves a limited number of statistical moments. The model parameters are sampled to satisfy these moments and collected in ensembles, similar to how parameters are sampled to fulfill probability distributions in RS.
The constraints of satisfying statistical moments constitute an infinite system of equations for the samples . It can formally be viewed as sampling () of the joint pdf ,
The infinite number of equations requires an infinite number of samples. However, it is implicitly assumed that relatively few moments are known and significantly influence the result of interest. Only a few moments then needs to be accurately represented by . Typically, and are estimated when models are identified [4,7]. In addition, the range or another higher diagonal moment can generally be determined from underlying physical constraints like stability. Clearly, any sampling rule must generate a fixed number of samples and create them simultaneously. The samples are consequently strongly dependent. One obvious sampling method is to solve Eq. 18 numerically for a sufficiently large number of samples , as in Eq. 16. Due to the strong non-linearities, this is quite difficult for a large number of moments but may be feasible for a few moments.
5.2. The excitation matrix
The UKF (section 4.2) utilizes DS with conservation of all first and second statistical moments. The invariance in its formulation allows for any additional ‘half’ unitary transformation , . This results in another equally valid matrix , since . Since the transformation is allowed and influences the result, the result of applying the UKF is not unique. The matrix condenses this invariance and provides practical means to manipulate the UKF ensemble. A key feature of is the absence of constraints on . That makes it possible to stretch ‘horizontally’ (as long as ). That corresponds to adding samples (sigma-points). The improved transformation (section 3.2) can be applied by also combining with . The square root of the covariance matrix will then read instead of ,
The samples are here collected in columns of the ensemble matrix . The matrix is diagonalized,, with the unitary transformation . The normalization factoris included in the excitation matrix to satisfy the correct covariance, , just as the factor was included in Eq. 17 (the ensemble is here expanded from to samples). The excitation matrix controls the sampling beyond the first and second moments, e.g. the range of the samples. Row of the matrix can be interpreted as deterministic samples of the pdf , assigned to canonical parameters in RS, see section 3.2. All ensembles will be described with a unique excitation matrix.
The adopted transformation of distorts all higher moments than the second. This mixing effect is indicated by the index defined in Eq. 11. To diagonalize large matrices many efficient techniques have been developed. This should not cause any difficulties even for , especially since usually is either very sparse, or rank deficient for models with many parameters.
5.3. Elimination of singular values
The rationale for applying the reduction to be presented is that any model is derived from a limited set of experiments, resulting in a usually moderate rank of . If the number of parameters is large, it is thus often (nearly) rank deficient. The widely practiced singular value decomposition (SVD)  may then be used to reduce the excitation matrix and hence the number of samples. The most general form of SVD cannot be used here since it renders an asymmetric decomposition , where and . Different matrices allow for decomposition of an arbitrary matrix. For the symmetric matrix , a symmetric SVD can be found with the less general eigenvalue decomposition , according to the spectral theorem . As is positive definite, all its eigenvalues fulfill the requirement of being positive. This is required to directly obtain a real-valued matrix square root .
The ensemble may now be reduced by elimination of singular values (ESV). Choose a threshold and remove row and column from and row from for all such that,
Proceeding as in many applications of SVD, this reduction (indicated by tilde below) will not change the result significantly, if is small enough. Accordingly, samples are eliminated using the alternative decomposition of the square root of ,
Unfortunately, the less distorting transformation of advocated in section 3.2 do not allow for rows of the matrix . The increase in distortion of indicated by is less important for the intended use though. Signal processing models with large numbers of parameters are typically non-parametric and usually describe samples of signals like impulse responses, or noise signals. The required LR property of the system then implies LP. The propagation of covariance is then linear and only the undistorted first and second moments need to be encoded.
5.4. Correlated sampling of non-parametric models
A major difference between parametric and non-parametric models is the dimensionality. A conceptual dissimilarity is that non-parametric models usually refer to correlated signals, rather than abstract model structures. The parameters may describe discrete samples of input noise , or an impulse response . A common parametric pole-zero model may contain 20 parameters, while a non-parametric model can be expressed in perhaps 1000 parameters. The ensembles of non-parametric models often need to be reduced drastically.
Due to limited resolution, the correlation times of any signal or impulse response is finite. Their ‘memory’ is thus finite so sample variations may be regenerated or repeated, as long as the time between repetitions exceeds the correlation time. This correlated sampling (CRS) provides efficient and accurate reduction of the ensembles. The minimal number of parameters is then set by the correlation time of the model. Most importantly, the size of the ensemble becomes independent of the size of the model (the length of the signal).
A finite correlation length of any model is normally inferred from the decay of its autocorrelation function, where denotes the lag and refers to a non-stationary variation. Here, a global will be defined through its -norm and determined for a relative truncation threshold (returns the minimizing argument),
If the model is expressed as a convolution of an impulse response and time-dependent white noise as in section 6.2,
By padding the model to an integer multiple of samples, it is always possible to choose an excitation matrix partitioned to block-diagonal form,
where is any allowed deterministic sub-ensemble. The factor accounts for the change from samples of to the samples of . By violating the normalization constraint , the size of the ensemble can be ‘compressed’ from to samples by moving all sub-matrices to the first block-column and skipping all zeros. The introduced constant drops out as ,
The consequence of violating the normalization constraint is that only a limited diagonal band of is correctly reproduced. If a non-parametric model of a signal is propagated through a system model with impulse response of correlation length this will nevertheless not result in any error of , as it is independent of all faulty elements , . To correctly evaluate though, the size of the sub-ensembles of correlated sampling must fulfill the stronger size constraint . The symmetry of convolutions implies a corresponding result when the non-parametric model describes the impulse response of the system, rather than a signal.
5.5. Combining covariance
A signal processing model generally includes both parametric and non-parametric sources of uncertainty. For instance, a device (parametric system model) may be fed with a signal corrupted with noise (non-parametric noise model). The question then arises how the two sources of uncertainty can be combined. For propagation of uncertainty through LP models, the combined covariance is given by the Gauss approximation formula ,
where is the propagated covariance of . This will seize to apply for non-LP models. There exists no general non-linear summation rule for propagated covariance. A method of summation can be given though, if different ensembles are combined as in RS.
To combine ensembles of parametric and non-parametric models , collect all parameters, , and diagonalize the enlarged covariance matrix, . Build with two blocks and use CRS (section 5.4) for the non-parametric model,
The scaling may cause a similar scaling problem as the factor in the UKF (section 4.2). Using extended excitation matrices these factors can be eliminated,
A disadvantage of this summation is that the same type of ensemble must be used for all parameters. Both alternatives combine the statistics of the two models non-linearly. The uncertainties are propagated and combined by evaluating the model for all samples and calculating the desired statistics, just as if the combined ensemble described one model.
5.6. Selected ensembles
The standard (STD) ensemble employed in the UKF (as defined in section 4.2) utilizes the perhaps simplest possible excitation matrix,
While the ultimate simplicity is its main advantage, the long maximal(!) range is its main disadvantage.
How far the reduction of samples might be driven is illustrated by the minimal simplex (SPX) ensemble,
where the operator performs classical Gram-Schmidt orthogonalization  and normalization of rows. The ensemble is constructed from half the STD ensemble, complemented with one sample to cancel the first moments. Since that violates the orthogonality of the rows of , must be applied to satisfy . The high efficiency of the SPX ensemble is tarnished by its large skewness, or . This may give considerable bias of propagated covariance for non-LP models, but is irrelevant for LP models.
The binary (BIN) ensemble has minimal range to guarantee allowed samples. By varying all parameters with an equal magnitude of one standard deviation in all samples, the diverging factor of the STD is eliminated. Its excitation matrix is fundamentally constructed from a standard binary array, with the difference that the allowed levels are instead of (see rows 1-3 in Eq. 32). It is then complemented with supplementary rows obtained in two ways, by cyclic shifting and mirror imaging,
Cyclic shifts are applied to all original rows except the first, by a quarter of their periodicity. Mirror imaging of a row is defined to change the sign of its second half and is applied to all original rows except the last two, and all shifted rows except the last. For instance, in Eq. 32 row 4 and 5 are shifted versions of row 2 and 3, while rows 6 and 7 are the mirror images of rows 1 and the shifted row 4. The supplementary rows reduce the size of the ensemble drastically with a corresponding improvement of the efficiency. For parameters, the size drops from roughly to samples. That size is acceptable in perspective of the samples of the most efficient SPX. Eventually though, the number of samples will grow too large. The BIN can thus only be applied to moderately sized models.
By no means, this brief survey exhausts all possible ensembles. Many criteria for selecting the most appropriate ensemble can be formulated. Here, the first and second moments, parameter ranges and efficiency were in focus.
6. Application — Modeling uncertainty of a dynamic device
The task is to simulate the response of an electrical device such as an amplifier or oscilloscope, in the presence of non-stationary correlated noise on its input. An uncertain LR CT model of the device and its parametric covariance is usually found by applying system identification techniques  on calibration measurements . Such a model of the system can be sampled into a digital filter and be described in the pole-zero form in Eq. 2. These standard steps will here be omitted. The system model will instead be assigned to a digital low-pass Butterworth filter, of order and cross-over frequency , being the Nyquist frequency and described by parameters . The complete correlations of complex-conjugated pole and zero pairs are eliminated by a transformation from to , giving system model parameters,
To be most general, the non-parametric input noise model is chosen to be correlated/colored and non-stationary. The noise parameter represents the noise level at time sample . Its generating signal  is a Dirac delta function , centered at time . The response of a system with impulse response will be . In matrix notation, , where . Hence,
Since the response is linear in noise parameters, it is sufficient to only capture .
6.1. Reference ensembles
Traditionally, any method for uncertainty propagation is evaluated by comparisons with the default method of linearization [10,16], and brute-force random sampling (MC)  as state-of-the-art. There are several drawbacks of this approach. Linearization is a coarse approximation for LP models and MC suffer from the difficulty of modeling dependencies and low efficiency. An alternative is to construct finite reference ensembles (REF) and by definition let them describe the truth. Their primary advantage is that the finite size of the REF makes it possible to propagate the uncertainty exactly, using all REF samples. A more or less arbitrary REF may be generated randomly, like any MC ensemble. All requirements are also automatically fulfilled since the REF is built of possible realizations. Also, the REF closes the loop as it makes it possible to compare ‘true’ and approximate samples directly on an equal footing (see Fig. 8). Even though the samples differ substantially, the resulting modeling uncertainties can be similar.
A plausible REF for the system model realized as a digital filter is created by randomly generating samples of parameters from uniform distributions, with listed in Fig. 3, top left. The joint pdf will have compact support , as required to guarantee stability. The mean is subtracted from all samples to remove the bias of the finite random ensemble, . The covariance of the REF will have a desirable more or less random variation for small values of . If the REF samples are arranged in columns of a matrix , (as ) . For , the strong correlations will expose the methods to severe tests with significant transformations . The mixing (Eq. 13) using transformation was considerable, but less than for , see caption Fig. 3. For the chosen REF, the resulting variations of poles and zeros are displayed in Fig. 3. The obtained variation of the parameters defined in Eq. 33 can be quantified with an averaged correlation index and standard deviation (Fig. 3, bottom left),
A REF signal for non-stationary correlated noise may conveniently be generated from an autoregressive process (AR) acting on time-dependent zero mean white noise , , where is the matrix of translated impulse responses for the AR process defined by parameters . Assigning a square wave time-dependence,
An explicit realization of the REF for the noise model is hence not needed. Specifically, a second order system with time parameters was chosen. The impulse responses of the AR noise system and the system model, and the variation of the noise model are illustrated in Fig. 3, bottom right.
The ‘true’ result given by the response for the REFs for the different test signals is shown in Fig. 4. The propagated noise variation differs substantially from the input square wave (top left) and is almost opposite in phase, due to the response time of about , see delay of (top, right and bottom). The signal distortion is strongly dependent on the input signal and decreases with increased regularity / differentiability. The propagated covariance has a more complex variation (top, right and bottom), as it is larger for the more regular Gaussian (bottom, right) than for the triangular pulse (bottom, left).
6.2. Deterministic sampling
The error of the scent and the standard deviation for the STD, SPX and BIN ensembles of the system model is displayed in Fig. 5, for all test signals. The low scent of the REF (left: thin, dotted) suggests the model is close to LP. Despite the relative errors are large they are quite small on an absolute scale. The SPX has the largest errors, for the scent as well as the variance. That is likely caused by its skewness being much larger than that of the REF. The BIN has the lowest errors and is thus the best approximate representation of the REF.
The errors might appear large, considering all ensembles are ‘correct’, i.e. correctly represent (typically) available accurate information (mean and covariance of parameters). The errors reflect ambiguities caused by the ubiquitous lack of information in signal processing, rather than inadequacies of DS. RS can only produce better results by making further assumptions.
The result of applying the ESV and the CRS methods to reduce the SPX ensemble for propagating the noise is displayed in Fig. 6. By choosing sufficiently low thresholds for elimination of singular values (ESV) and for truncation of the correlation lengths (see Eqs. 20,22), the errors can be made arbitrarily low. As the reduction will decrease accordingly, there is a trade-off between accuracy and efficiency. For the chosen values, CRS is about twice as accurate and twice as efficient as ESV. In contrast to ESV, the number of samples for the CRS method is independent of the number of noise samples. The computational cost thus increases linearly with the length of the noise signal for CRS but quadratically (approximately) for ESV. For ESV to be most efficient, the model covariance needs to be strongly rank deficient. That is not as unlikely as it might appear, since the model usually is derived from a limited amount of experimental results.
The summation of the noise and the model covariance is illustrated in Fig. 7. The propagation of the covariance of the system model is not LP. The quadratic summation rule (Eq. 27), or Gauss approximation formula , is therefore not applicable. Nevertheless, the low scent (Fig. 5, left) suggests that both propagations are close to LP. The summation error is hence finite, but quite small. It differs qualitatively from both contributions, indicating that the summation is non-trivial.
Finally, the samples of one pole of the derived ensembles are compared to the reference samples of the REF in Fig. 8. The limit of stability is included to illustrate how close the samples are to be physically forbidden. The construction of the different ensembles is apparent, even though the transformation distorts the scatter plots (sections 3.2, 5.2), and tilts the principal axes (lines). The samples of the REF are almost evenly distributed. Only four samples of the STD, labelled , deviate significantly from a dense central cluster, as described by the excitation matrix (Eq. 30). It also is evident that SPX originates from half the STD. A small translation required to achieve the correct mean is discernible, while the Gram-Schmidt orthogonalization renders a minor rotation and distortion. The BIN contains comparable variations in all samples and thus has no central cluster and its samples are repelled from the principal directions (lines). The statistical differences to the REF refer to the shape of the joint pdf. Choosing the best ensemble is thus equivalent of selecting the most appropriate pdf in RS. The BIN seems to resemble the REF scatter plot the most, as verified by its low errors in Fig. 5.
Deterministic sampling remains controversial  while random sampling has qualified as a preferred state-of-the-art method for propagating uncertainty. Both result in finite statistical  ensembles, which are approximate finite representations of the primary statistical models. Their sampling strategies and convergence rates are dramatically different. While deterministic sampling humbly aims at representing the most relevant and best known statistical information, random sampling targets complete control of all features of the ensemble. Such detailed information is rarely known and must instead be more or less blindly assigned. The inevitable consequence is that critical computational resources are spent on propagating, at best, vaguely known details. The numerical power of modern computers is better spent on refinements of the signal processing model (longer time series, higher sampling rates, larger systems etc.). Refined methods of random sampling have therefore been proposed which either simplifies the model, or improve the sampling distributions. Compared to deterministic sampling though, their convergence rates remain low.
It is easy to confuse deterministic sampling with experimental design and optimization . Even though any sample could be a possible outcome of an experiment, deterministic ensembles represent rather than realize (as random ensembles) statistical distributions. Instead of associating a joint distribution to the parameters of an uncertain model, it is possible to directly represent their statistics with a deterministic ensemble. That would eliminate the need of interpreting abstract distributions and result in complete reproducibility. The critical choice of ensemble would be assigned once and for all in the calibration experiment, with no further need of approximation.
The use of excitation matrices made it possible to construct universal generic ensembles. The efficiency of the minimal SPX ensemble is indeed high but so is also its third moment. While the STD maximizes the range of each parameter, the BIN minimizes it by varying all parameters in all samples. The STD is the simplest while the SPX is the most efficient ensemble. In the example, the BIN was most accurate. For non-parametric models with many parameters, reduction of samples may be required. Elimination of singular values (ESV) and correlated sampling (CRS) were two such techniques. The presented ensembles are not to be associated to random sampling as a method. They are nothing but a few examples of deterministic sampling, likely the best ensembles are yet to be discovered.
It is indeed challenging but also rewarding to find novel deterministic sampling strategies. Once the sampling rules are found, the application is just as simple as random sampling, but usually much more efficient. Deterministic sampling is one of very few methods capable of non-linear propagation of uncertainty through large signal processing models.