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
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
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
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
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
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
3.2. Brute force Monte Carlo
Monte Carlo (MC) methods [8-9], or
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
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
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
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
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
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
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
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
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
The samples are here collected in columns of the ensemble matrix . The matrix is diagonalized,, with the unitary transformation . The normalization factor is included in the
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
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
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
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
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
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
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.
Kay S. Fundamentals of Statistical signal processing: Estimation Theory. New Jersey: Prentice Hall; 1993.
Hessling JP. Propagation of dynamic measurement uncertainty. Meas. Sci. Technol. 2011; 22 (10) 105105 (13pp).
Hessling JP. Integration of digital filters and measurements. In: Márquez FPG. (ed.) Digital Filters. Rijeka: InTech; 2011. p123-154. Available from http://www.intechopen.com/books/digital-filters/integration-of-digital-filters-and-measurements (accessed 4 July 2012).
Pintelon R, Schoukens J. System Identification: A Frequency Domain Approach. Piscataway, New Jersey: IEEE Press; 2001.
Witteveen JAS. Efficient and Robust Uncertainty Quantification for Computational Fluid Dynamics and Fluid-Structure Interaction. PhD thesis. Delft University of Technology; 2009.
Hale PD, Dienstfrey A, Wang JCM, Williams DF, Lewandowski A, Keenan DA, Clement TS. Traceable Waveform Calibration With a Covariance-Based Uncertainty Analysis. IEEE Trans. Instrum. Meas. 2009; 58 (10) 3554-3568.
Hessling JP. Metrology for non-stationary dynamic measurements. In: Sharma MK. (ed.) Advances in Measurement systems. Vukovar: InTech; 2010. p. 221-256. Available from http://www.intechopen.com/books/advances-in-measurement-systems/metrology-for-non-stationary-dynamic-measurements (accessed 4 July 2012).
Metropolis N, Ulam S. The Monte Carlo Method. Journal of the American Statistical Association 1949; 44 (247) 335-341.
Rubenstein RY, Kroese DP. Simulation and the Monte Carlo Method, 2nd Ed. New York: John Wiley & Sons Inc.; 2007.
Hessling JP. A novel method of evaluating dynamic measurement uncertainty utilizing digital filters. Meas. Sci. Technol. 2009; 20 (5) 055106 (11pp).
Råde L, Westergren B. Beta Mathematics Handbook, 2nd Ed. Lund, Sweden: Studentlitteratur; 1990.
Hessling JP, Svensson T. Propagation of uncertainty by sampling on confidence boundaries, accepted for publication in International Journal for Uncertainty Quantification.
Hessling JP. Deterministic sampling for propagating model covariance, submitted for publication.
Lovett T. Polynomial Chaos Simulation of Analog and Mixed-Signal Systems: Theory, Modeling method, Application. Saarbrucken: Lambert Academic Publishing; 2010.
Li H, Zhang D. Probabilistic collocation method for flow in porous media: Comparisons with other stochastic methods. Water Resources Research 2007; 43 W09409 (13 pp).
ISO GUM. Guide to the Expression of Uncertainty in Measurement. Geneva: International Organisation for Standardisation; 1995.
Metropolis N. The Beginning of the Monte Carlo Method. Los Alamos Science special issue 1987; 15 125-130.
Helton J, Davis L. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliability Engineering and System Safety 2003; 81 23-69.
Björk Å. Numerical methods for least squares problems. Philadelphia: Siam; 1996.
Wikipedia: http://en.wikipedia.org/wiki/Latin_square (accessed 3 July 2012):
Julier S, Uhlmann J, Durrant-Whyte HF. A new approach for filtering nonlinear systems. Proc IEEE American Control Conference June 21-23 1995; 1628-1632.
Julier S, Uhlmann J. Unscented filtering and nonlinear estimation. Proceeding IEEE March 2004; 92 (3) 401-422.
Simon D. Optimal State Estimation: Kalman, and non-linear approaches. New Jersey: Wiley; 2006.
Matlab with Signal Processing Toolbox, The Mathworks, Inc.
Julier S, Uhlmann J. The scaled unscented transformation, Proceedings of the IEEE American Control Conference 8-10 May 2002; 4555-4559.
Hessling JP. Non-linear propagation and summation of covariance using deterministic sampling, in preparation.
Gustafsson F, Hendeby G. Some Relations between Extended and Unscented Kalman Filters. IEEE Trans. sign. proc. 2012; 60 (2) 545-555.
Fischer RA. Statistical Methods, Experimental design and Scientific Inference. New York: Oxford University Press; 1990.