Evaluation of uncertainties¶
The evaluation of uncertainties is a fundamental part of the measurement
analysis in metrology. The analysis of dynamic measurements typically
involves methods from signal processing, such as digital filtering or
application of the discrete Fourier transform (DFT). For most tasks, methods
are readily available, for instance, as part of scipy.signals
. This
module of PyDynamic provides the corresponding methods for the evaluation of
uncertainties.
Uncertainty evaluation for the DFT¶
The PyDynamic.uncertainty.propagate_DFT
module implements methods for
the propagation of uncertainties in the application of the DFT, inverse DFT,
deconvolution and multiplication in the frequency domain, transformation from
amplitude and phase to real and imaginary parts and vice versa.
 The corresponding scientific publications is
 S. Eichstädt und V. Wilkens GUM2DFT — a software tool for uncertainty evaluation of transient signals in the frequency domain. Measurement Science and Technology, 27(5), 055001, 2016. [DOI: 10.1088/09570233/27/5/055001]
This module contains the following functions:
 GUM_DFT: Calculation of the DFT of the time domain signal x and propagation of the squared uncertainty Ux associated with the time domain sequence x to the real and imaginary parts of the DFT of x
 GUM_iDFT: GUM propagation of the squared uncertainty UF associated with the DFT values F through the inverse DFT
 GUM_DFTfreq: Return the Discrete Fourier Transform sample frequencies
 DFT_transferfunction: Calculation of the transfer function H = Y/X in the frequency domain with X being the Fourier transform of the system’s input signal and Y that of the output signal
 DFT_deconv: Deconvolution in the frequency domain
 DFT_multiply: Multiplication in the frequency domain
 AmpPhase2DFT: Transformation from magnitude and phase to real and imaginary parts
 DFT2AmpPhase: Transformation from real and imaginary parts to magnitude and phase
 AmpPhase2Time: Transformation from amplitude and phase to time domain
 Time2AmpPhase: Transformation from time domain to amplitude and phase

PyDynamic.uncertainty.propagate_DFT.
GUM_DFT
(x, Ux, N=None, window=None, CxCos=None, CxSin=None, returnC=False, mask=None)[source]¶ Calculation of the DFT with propagation of uncertainty
Calculation of the DFT of the time domain signal x and propagation of the squared uncertainty Ux associated with the time domain sequence x to the real and imaginary parts of the DFT of x.
Parameters:  x (numpy.ndarray of shape (M,)) – vector of time domain signal values
 Ux (numpy.ndarray) – covariance matrix associated with x, shape (M,M) or vector of squared standard uncertainties, shape (M,) or noise variance as float
 N (int, optional) – length of time domain signal for DFT; N>=len(x)
 window (numpy.ndarray, optional of shape (M,)) – vector of the time domain window values
 CxCos (numpy.ndarray, optional) – cosine part of sensitivity matrix
 CxSin (numpy.ndarray, optional) – sine part of sensitivity matrix
 returnC (bool, optional) – if true, return sensitivity matrix blocks for later use
 mask (ndarray of dtype bool) – calculate DFT values and uncertainties only at those frequencies where mask is True
Returns:  F (numpy.ndarray) – vector of complex valued DFT values or of its real and imaginary parts
 UF (numpy.ndarray) – covariance matrix associated with real and imaginary part of F
References
 Eichstädt and Wilkens [Eichst2016]

PyDynamic.uncertainty.propagate_DFT.
GUM_iDFT
(F, UF, Nx=None, Cc=None, Cs=None, returnC=False)[source]¶ GUM propagation of the squared uncertainty UF associated with the DFT values F through the inverse DFT
The matrix UF is assumed to be for real and imaginary part with blocks: UF = [[u(R,R), u(R,I)],[u(I,R),u(I,I)]] and real and imaginary part obtained from calling rfft (DFT for realvalued signal)
Parameters:  F (np.ndarray of shape (2M,)) – vector of real and imaginary parts of a DFT result
 UF (np.ndarray of shape (2M,2M)) – covariance matrix associated with real and imaginary parts of F
 Nx (int, optional) – number of samples of iDFT result
 Cc (np.ndarray, optional) – cosine part of sensitivities (without scaling factor 1/N)
 Cs (np.ndarray, optional) – sine part of sensitivities (without scaling factor 1/N)
 returnC (if true, return sensitivity matrix blocks (without scaling factor 1/N)) –
Returns:  x (np.ndarray) – vector of time domain signal values
 Ux (np.ndarray) – covariance matrix associated with x
References
 Eichstädt and Wilkens [Eichst2016]

PyDynamic.uncertainty.propagate_DFT.
GUM_DFTfreq
(N, dt=1)[source]¶ Return the Discrete Fourier Transform sample frequencies
Parameters:  N (int) – window length
 dt (float) – sample spacing (inverse of sampling rate)
Returns: f – Array of length
n//2 + 1
containing the sample frequenciesReturn type: ndarray

PyDynamic.uncertainty.propagate_DFT.
DFT_transferfunction
(X, Y, UX, UY)[source]¶ Calculation of the transfer function H = Y/X in the frequency domain
Calculate the transfer function with X being the Fourier transform of the system’s input signal and Y that of the output signal.
Parameters:  X (np.ndarray) – real and imaginary parts of the system’s input signal
 Y (np.ndarray) – real and imaginary parts of the system’s output signal
 UX (np.ndarray) – covariance matrix associated with X
 UY (np.ndarray) – covariance matrix associated with Y
Returns:  H (np.ndarray) – real and imaginary parts of the system’s frequency response
 UH (np.ndarray) – covariance matrix associated with H
This function only calls DFT_deconv.

PyDynamic.uncertainty.propagate_DFT.
DFT_deconv
(H, Y, UH, UY)[source]¶ Deconvolution in the frequency domain
GUM propagation of uncertainties for the deconvolution X = Y/H with Y and H being the Fourier transform of the measured signal and of the system’s impulse response, respectively. This function returns the covariance matrix as a tuple of blocks if too large for complete storage in memory.
Parameters:  H (np.ndarray of shape (2M,)) – real and imaginary parts of frequency response values (N an even integer)
 Y (np.ndarray of shape (2M,)) – real and imaginary parts of DFT values
 UH (np.ndarray of shape (2M,2M)) – covariance matrix associated with H
 UY (np.ndarray of shape (2M,2M)) – covariance matrix associated with Y
Returns:  X (np.ndarray of shape (2M,)) – real and imaginary parts of DFT values of deconv result
 UX (np.ndarray of shape (2M,2M)) – covariance matrix associated with real and imaginary part of X
References
 Eichstädt and Wilkens [Eichst2016]

PyDynamic.uncertainty.propagate_DFT.
DFT_multiply
(Y, F, UY, UF=None)[source]¶ Multiplication in the frequency domain
GUM uncertainty propagation for multiplication in the frequency domain, where the second factor F may have an associated uncertainty. This method can be used, for instance, for the application of a lowpass filter in the frequency domain or the application of deconvolution as a multiplication with an inverse of known uncertainty.
Parameters:  Y (np.ndarray of shape (2M,)) – real and imaginary parts of the first factor
 F (np.ndarray of shape (2M,)) – real and imaginary parts of the second factor
 UY (np.ndarray either shape (2M,) or shape (2M,2M)) – covariance matrix or squared uncertainty associated with Y
 UF (np.ndarray of shape (2M,2M)) – covariance matrix associated with F (optional), default is None
Returns:  YF (np.ndarray of shape (2M,)) – the product of Y and F
 UYF (np.ndarray of shape (2M,2M)) – the uncertainty associated with YF

PyDynamic.uncertainty.propagate_DFT.
AmpPhase2DFT
(A, P, UAP, keep_sparse=False)[source]¶ Transformation from magnitude and phase to real and imaginary parts
Calculate the vector F=[real,imag] and propagate the covariance matrix UAP associated with [A, P]
Parameters:  A (np.ndarray of shape (N,)) – vector of magnitude values
 P (np.ndarray of shape (N,)) – vector of phase values (in radians)
 UAP (np.ndarray of shape (2N,2N)) – covariance matrix associated with (A,P) or vector of squared standard uncertainties [u^2(A),u^2(P)]
 keep_sparse (bool, optional) – whether to transform sparse matrix to numpy array or not
Returns:  F (np.ndarray) – vector of real and imaginary parts of DFT result
 UF (np.ndarray) – covariance matrix associated with F

PyDynamic.uncertainty.propagate_DFT.
DFT2AmpPhase
(F, UF, keep_sparse=False, tol=1.0, return_type='separate')[source]¶ Transformation from real and imaginary parts to magnitude and phase
Calculate the matrix U_AP = [[U1,U2],[U2^T,U3]] associated with magnitude and phase of the vector F=[real,imag] with associated covariance matrix U_F=[[URR,URI],[URI^T,UII]]
Parameters:  F (np.ndarray of shape (2M,)) – vector of real and imaginary parts of a DFT result
 UF (np.ndarray of shape (2M,2M)) – covariance matrix associated with F
 keep_sparse (bool, optional) – if true then UAP will be sparse if UF is onedimensional
 tol (float, optional) – lower bound for A/uF below which a warning will be issued concerning unreliable results
 return_type (str, optional) – If “separate” then magnitude and phase are returned as separate arrays. Otherwise the array [A, P] is returned
If return_type is separate:
Returns:  A (np.ndarray) – vector of magnitude values
 P (np.ndarray) – vector of phase values in radians, in the range [pi, pi]
 UAP (np.ndarray) – covariance matrix associated with (A,P)
Otherwise:
Returns:  AP (np.ndarray) – vector of magnitude and phase values
 UAP (np.ndarray) – covariance matrix associated with AP

PyDynamic.uncertainty.propagate_DFT.
AmpPhase2Time
(A, P, UAP)[source]¶ Transformation from amplitude and phase to time domain
GUM propagation of covariance matrix UAP associated with DFT amplitude A and phase P to the result of the inverse DFT. Uncertainty UAP is assumed to be given for amplitude and phase with blocks: UAP = [[u(A,A), u(A,P)],[u(P,A),u(P,P)]]
Parameters:  A (np.ndarray of shape (N,)) – vector of amplitude values
 P (np.ndarray of shape (N,)) – vector of phase values (in rad)
 UAP (np.ndarray of shape (2N,2N)) – covariance matrix associated with [A,P]
Returns:  x (np.ndarray) – vector of time domain values
 Ux (np.ndarray) – covariance matrix associated with x

PyDynamic.uncertainty.propagate_DFT.
Time2AmpPhase
(x, Ux)[source]¶ Transformation from time domain to amplitude and phase
Parameters:  x (np.ndarray of shape (N,)) – time domain signal
 Ux (np.ndarray of shape (N,N)) – squared uncertainty associated with x
Returns:  A (np.ndarray) – amplitude values
 P (np.ndarray) – phase values
 UAP (np.ndarray) – covariance matrix associated with [A,P]

PyDynamic.uncertainty.propagate_DFT.
Time2AmpPhase_multi
(x, Ux, selector=None)[source]¶ Transformation from time domain to amplitude and phase
Perform transformation for a set of M signals of the same type.
Parameters:  x (np.ndarray of shape (M,N)) – M time domain signals of length N
 Ux (np.ndarray of shape (M,)) – squared standard deviations representing noise variances of the signals x
 selector (np.ndarray of shape (L,), optional) – indices of amplitude and phase values that should be returned; default is 0:N1
Returns:  A (np.ndarray of shape (M,N)) – amplitude values
 P (np.ndarray of shape (M,N)) – phase values
 UAP (np.ndarray of shape (M, 3N)) – diagonals of the covariance matrices: [diag(UPP), diag(UAA), diag(UPA)]
Uncertainty evaluation for digital filtering¶
This module contains functions for the propagation of uncertainties through the application of a digital filter using the GUM approach.
This modules contains the following functions:
 FIRuncFilter: Uncertainty propagation for signal y and uncertain FIR filter theta
 IIRuncFilter: Uncertainty propagation for the signal x and the uncertain IIR filter (b,a)
# Note: The ElsterLink paper for FIR filters assumes that the autocovariance is known and that noise is stationary!

PyDynamic.uncertainty.propagate_filter.
FIRuncFilter
(y, sigma_noise, theta, Utheta=None, shift=0, blow=None, kind='corr')[source]¶ Uncertainty propagation for signal y and uncertain FIR filter theta
Parameters:  y (np.ndarray) – filter input signal
 sigma_noise (float or np.ndarray) – float: standard deviation of white noise in y 1Darray: interpretation depends on kind
 theta (np.ndarray) – FIR filter coefficients
 Utheta (np.ndarray) – covariance matrix associated with theta
 shift (int) – time delay of filter output signal (in samples)
 blow (np.ndarray) – optional FIR lowpass filter
 kind (string) – only meaningfull in combination with isinstance(sigma_noise, numpy.ndarray) “diag”: pointwise standard uncertainties of nonstationary white noise “corr”: single sided autocovariance of stationary (colored/corrlated) noise (default)
Returns:  x (np.ndarray) – FIR filter output signal
 ux (np.ndarray) – pointwise uncertainties associated with x
References
 Elster and Link 2008 [Elster2008]
See also

PyDynamic.uncertainty.propagate_filter.
IIRuncFilter
(x, noise, b, a, Uab)[source]¶ Uncertainty propagation for the signal x and the uncertain IIR filter (b,a)
Parameters:  x (np.ndarray) – filter input signal
 noise (float) – signal noise standard deviation
 b (np.ndarray) – filter numerator coefficients
 a (np.ndarray) – filter denominator coefficients
 Uab (np.ndarray) – covariance matrix for (a[1:],b)
Returns:  y (np.ndarray) – filter output signal
 Uy (np.ndarray) – uncertainty associated with y
References
 Link and Elster [Link2009]
Monte Carlo methods for digital filtering¶
The propagation of uncertainties via the FIR and IIR formulae alone does not enable the derivation of credible intervals, because the underlying distribution remains unknown. The GUMS2 Monte Carlo method provides a reference method for the calculation of uncertainties for such cases.
This module implements Monte Carlo methods for the propagation of uncertainties for digital filtering.
This module contains the following functions:
 MC: Standard Monte Carlo method for application of digital filter
 SMC: Sequential Monte Carlo method with reduced computer memory requirements
 UMC: Update Monte Carlo method for application of digital filters with reduced computer memory requirements
 UMC_generic: Update Monte Carlo method with reduced computer memory requirements

PyDynamic.uncertainty.propagate_MonteCarlo.
MC
(x, Ux, b, a, Uab, runs=1000, blow=None, alow=None, return_samples=False, shift=0, verbose=True)[source]¶ Standard Monte Carlo method
Monte Carlo based propagation of uncertainties for a digital filter (b,a) with uncertainty matrix \(U_{\theta}\) for \(\theta=(a_1,\ldots,a_{N_a},b_0,\ldots,b_{N_b})^T\)
Parameters:  x (np.ndarray) – filter input signal
 Ux (float or np.ndarray) – standard deviation of signal noise (float), pointwise standard uncertainties or covariance matrix associated with x
 b (np.ndarray) – filter numerator coefficients
 a (np.ndarray) – filter denominator coefficients
 Uab (np.ndarray) – uncertainty matrix \(U_\theta\)
 runs (int,optional) – number of Monte Carlo runs
 return_samples (bool, optional) – whether samples or mean and std are returned
If
return_samples
isFalse
, the method returns:Returns:  y (np.ndarray) – filter output signal
 Uy (np.ndarray) – uncertainty associated with
Otherwise the method returns:
Returns: Y – array of Monte Carlo results Return type: np.ndarray References
 Eichstädt, Link, Harris and Elster [Eichst2012]

PyDynamic.uncertainty.propagate_MonteCarlo.
SMC
(x, noise_std, b, a, Uab=None, runs=1000, Perc=None, blow=None, alow=None, shift=0, return_samples=False, phi=None, theta=None, Delta=0.0)[source]¶ Sequential Monte Carlo method
Sequential Monte Carlo propagation for a digital filter (b,a) with uncertainty matrix \(U_{\theta}\) for \(\theta=(a_1,\ldots,a_{N_a},b_0,\ldots,b_{N_b})^T\)
Parameters:  x (np.ndarray) – filter input signal
 noise_std (float) – standard deviation of signal noise
 b (np.ndarray) – filter numerator coefficients
 a (np.ndarray) – filter denominator coefficients
 Uab (np.ndarray) – uncertainty matrix \(U_\theta\)
 runs (int, optional) – number of Monte Carlo runs
 Perc (list, optional) – list of percentiles for quantile calculation
 blow (np.ndarray) – optional lowpass filter numerator coefficients
 alow (np.ndarray) – optional lowpass filter denominator coefficients
 shift (int) – integer for time delay of output signals
 return_samples (bool, otpional) – whether to return y and Uy or the matrix Y of MC results
 theta (phi,) – parameters for AR(MA) noise model \(\epsilon(n) = \sum_k \phi_k\epsilon(nk) + \sum_k \theta_k w(nk) + w(n)\) with \(w(n)\sim N(0,noise_std^2)\)
 Delta (float,optional) – upper bound on systematic error of the filter
If
return_samples
isFalse
, the method returns:Returns:  y (np.ndarray) – filter output signal (Monte Carlo mean)
 Uy (np.ndarray) – uncertainties associated with y (Monte Carlo pointwise std)
 Quant (np.ndarray) – quantiles corresponding to percentiles
Perc
(if notNone
)
Otherwise the method returns:
Returns: Y – array of all Monte Carlo results Return type: np.ndarray References
 Eichstädt, Link, Harris, Elster [Eichst2012]

PyDynamic.uncertainty.propagate_MonteCarlo.
UMC
(x, b, a, Uab, runs=1000, blocksize=8, blow=1.0, alow=1.0, phi=0.0, theta=0.0, sigma=1, Delta=0.0, runs_init=100, nbins=1000, credible_interval=0.95)[source]¶ Batch Monte Carlo for filtering using update formulae for mean, variance and (approximated) histogram. This is a wrapper for the UMC_generic function, specialised on filters
Parameters:  x (np.ndarray, shape (nx, )) – filter input signal
 b (np.ndarray, shape (nbb, )) – filter numerator coefficients
 a (np.ndarray, shape (naa, )) – filter denominator coefficients, normalization (a[0] == 1.0) is assumed
 Uab (np.ndarray, shape (naa + nbb  1, )) – uncertainty matrix \(U_\theta\)
 runs (int, optional) – number of Monte Carlo runs
 blocksize (int, optional) – how many samples should be evaluated for at a time
 blow (float or np.ndarray, optional) – filter coefficients of optional low pass filter
 alow (float or np.ndarray, optional) – filter coefficients of optional low pass filter
 phi (np.ndarray, optional,) – see misc.noise.ARMA noise model
 theta (np.ndarray, optional) – see misc.noise.ARMA noise model
 sigma (float, optional) – see misc.noise.ARMA noise model
 Delta (float, optional) – upper bound of systematic correction due to regularisation (assume uniform distribution)
 runs_init (int, optional) – how many samples to evaluate to form initial guess about limits
 nbins (int, list of int, optional) – number of bins for histogram
 credible_interval (float, optional) – must be in [0,1] central credible interval size
By default, phi, theta, sigma are chosen such, that N(0,1)noise is added to the input signal.
Returns:  y (np.ndarray) – filter output signal
 Uy (np.ndarray) – uncertainty associated with
 y_cred_low (np.ndarray) – lower boundary of credible interval
 y_cred_high (np.ndarray) – upper boundary of credible interval
 happr (dict) – dictionary keys: given nbin dictionary values: binedges val[“binedges”], bincounts val[“bincounts”]
References
 Eichstädt, Link, Harris, Elster [Eichst2012]
 ported to python in 201908 from matlabversion of Sascha Eichstaedt (PTB) from 20111012
 copyright on updating formulae parts is by Peter Harris (NPL)

PyDynamic.uncertainty.propagate_MonteCarlo.
UMC_generic
(draw_samples, evaluate, runs=100, blocksize=8, runs_init=10, nbins=100, return_samples=False, n_cpu=1)[source]¶ Generic Batch Monte Carlo using update formulae for mean, variance and (approximated) histogram. Assumes that the input and output of evaluate are numeric vectors (but not necessarily of same dimension). If the output of evaluate is multidimensional, it will be flattened into 1D.
Parameters:  draw_samples (function(int nDraws)) – function that draws nDraws from a given distribution / population needs to return a list of (multi dimensional) numpy.ndarrays
 evaluate (function(sample)) – function that evaluates a sample and returns the result needs to return a (multi dimensional) numpy.ndarray
 runs (int, optional) – number of Monte Carlo runs
 blocksize (int, optional) – how many samples should be evaluated for at a time
 runs_init (int, optional) – how many samples to evaluate to form initial guess about limits
 nbins (int, list of int, optional) – number of bins for histogram
 return_samples (bool, optional) – see returnvalue of documentation
 n_cpu (int, optional) – number of CPUs to use for multiprocessing, defaults to all available CPUs
Example
draw samples from multivariate normal distribution:
draw_samples = lambda size: np.random.multivariate_normal(x, Ux, size)
build a function, that only accepts one argument by masking additional kwargs:
evaluate = functools.partial(_UMCevaluate, nbb=b.size, x=x, Delta=Delta, phi=phi, theta=theta, sigma=sigma, blow=blow, alow=alow)
evaluate = functools.partial(bigFunction, **dict_of_kwargs)
By default the method
Returns:  y (np.ndarray) – mean of flattened/raveled simulation output i.e.: y = np.ravel(evaluate(sample))
 Uy (np.ndarray) – covariance associated with y
 happr (dict) – dictionary of binedges and bincounts
 output_shape (tuple) – shape of the unraveled simulation output can be used to reshape y and np.diag(Uy) into original shape
If
return_samples
isTrue
, the method additionally returns all evaluated samples. This should only be done for testing and debugging reasons, as this removes all memoryimprovements of the UMCmethod.Returns: sims – dict of samples and corresponding results of every evaluated simulation samples and results are saved in their original shape Return type: dict References
 Eichstädt, Link, Harris, Elster [Eichst2012]