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, the discrete
Fourier transform (DFT), or simple tasks like interpolation. For most of these tasks,
methods are readily available, for instance, as part of scipy.signal
. This
module of PyDynamic provides the corresponding methods for the evaluation of
uncertainties.
The package consists of the following modules:
PyDynamic.uncertainty.propagate_DFT
: Uncertainty evaluation for the DFTPyDynamic.uncertainty.propagate_DWT
: Uncertainty evaluation for the DWTPyDynamic.uncertainty.propagate_convolution
: Uncertainty evaluation for convolutionsPyDynamic.uncertainty.propagate_filter
: Uncertainty evaluation for digital filteringPyDynamic.uncertainty.propagate_MonteCarlo
: Monte Carlo methods for digital filteringPyDynamic.uncertainty.interpolate
: Uncertainty evaluation for interpolation
Uncertainty evaluation for convolutions
This module assists in uncertainty propagation for the convolution operation
The convolution operation is a common operation in signal and data processing. Convolving signals is mathematically similar to a filter application.
This module contains the following function:
convolve_unc()
: Convolution with uncertainty propagation based on FIR-filter
- PyDynamic.uncertainty.propagate_convolution.convolve_unc(x1, U1, x2, U2, mode='full')[source]
Discrete convolution of two signals with uncertainty propagation
This function supports the convolution modes of
numpy.convolve()
andscipy.ndimage.convolve1d()
.Note
The option to provide the uncertainties as 1D-arrays of standard uncertainties is given for convenience only. It does not result in any performance benefits, as they are internally just converted into a diagonal covariance matrix. Moreover, the output will always be a full covariance matrix (and will almost always have off-diagonal entries in practical scenarios).
- Parameters:
x1 (np.ndarray, (N,)) – first input signal
U1 (np.ndarray, (N,) or (N, N)) –
1D-array: standard uncertainties associated with x1 (corresponding to uncorrelated entries of x1)
2D-array: full 2D-covariance matrix associated with x1
None: corresponds to a fully certain signal x1, results in more efficient calculation (compared to using np.zeros(…))
x2 (np.ndarray, (M,)) – second input signal
U2 (np.ndarray, (M,) or (M, M)) –
1D-array: standard uncertainties associated with x2 (corresponding to uncorrelated entries of x2)
2D-array: full 2D-covariance matrix associated with x2
None: corresponds to a fully certain signal x2, results in more efficient calculation (compared to using np.zeros(…))
mode (str, optional) –
numpy.convolve()
-modes:full: len(y) == N+M-1 (default)
valid: len(y) == max(M, N) - min(M, N) + 1
same: len(y) == max(M, N) (value+covariance are padded with zeros)
scipy.ndimage.convolve1d()
-modes:nearest: len(y) == N (value+covariance are padded with by stationary assumption)
reflect: len(y) == N
mirror: len(y) == N
- Returns:
conv (np.ndarray) – convoluted output signal
Uconv (np.ndarray) – full 2D-covariance matrix of y
References
Uncertainty evaluation for the DFT
Functions for the propagation of uncertainties in the application of the DFT
The PyDynamic.uncertainty.propagate_DFT
module implements functions 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/0957-0233/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 xGUM_iDFT()
: GUM propagation of the squared uncertainty UF associated with the DFT values F through the inverse DFTGUM_DFTfreq()
: Return the Discrete Fourier Transform sample frequenciesDFT_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 signalDFT_deconv()
: Deconvolution in the frequency domainDFT_multiply()
: Multiplication in the frequency domainAmpPhase2DFT()
: Transformation from magnitude and phase to real and imaginary partsDFT2AmpPhase()
: Transformation from real and imaginary parts to magnitude and phaseAmpPhase2Time()
: Transformation from amplitude and phase to time domainTime2AmpPhase()
: Transformation from time domain to amplitude and phase
- PyDynamic.uncertainty.propagate_DFT.AmpPhase2DFT(A: ndarray, P: ndarray, UAP: ndarray, keep_sparse: bool = False) Tuple[ndarray, ndarray] [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) or of shape (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 of shape (2N,)) – vector of real and imaginary parts of DFT result
UF (np.ndarray of shape (2N,2N)) – covariance matrix associated with F
- Raises:
ValueError – If dimensions of A, P and UAP do not match.
- PyDynamic.uncertainty.propagate_DFT.AmpPhase2Time(A: ndarray, P: ndarray, UAP: ndarray) Tuple[ndarray, ndarray] [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 of shape (N, )) – vector of time domain values
Ux (np.ndarray of shape (2N, 2N)) – covariance matrix associated with x
- Raises:
ValueError – If dimension of UAP is not even.
- PyDynamic.uncertainty.propagate_DFT.DFT2AmpPhase(F: ndarray, UF: ndarray, keep_sparse: bool = False, tol: float = 1.0, return_type: str = 'separate') Tuple[ndarray, ndarray] | Tuple[ndarray, ndarray, ndarray] [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 one-dimensional
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 A and P. Otherwise the list [A, P] is returned
separate (If return_type ==) –
- Returns:
A (np.ndarray) – vector of magnitude values
P (np.ndarray) – vector of phase values in radians, in the range [-pi, pi], but only present if
return_type = 'separate'
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.DFT_deconv(H: ndarray, Y: ndarray, UH: ndarray, UY: ndarray) Tuple[ndarray, Tuple[ndarray, ndarray, ndarray] | ndarray] [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 (M an even integer)
Y (np.ndarray of shape (2M,)) – real and imaginary parts of DFT values
UH (np.ndarray of shape (2M,2M) or (2M,)) – full covariance or diagonal of the covariance matrix associated with H
UY (np.ndarray of shape (2M,2M) or (2M,)) – full covariance or diagonal of the 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) or 3-tuple of np.ndarray of shape (M,M)) – Covariance matrix associated with real and imaginary part of X. If the matrix fully assembled does not fit the memory, we return the auto-covariance for the real parts
URRX
and the imaginary partsUIIX
and the covariance between the real and imaginary partsURIX
as separatenp.ndarrays
arranged as follows:(URRX, URIX, UIIX)
References
Eichstädt and Wilkens [Eichst2016]
- Raises:
ValueError – If dimensions of H, Y, UY and UH do not match accordingly.
- PyDynamic.uncertainty.propagate_DFT.DFT_multiply(Y: ndarray, F: ndarray, UY: ndarray, UF: ndarray | None = None) Tuple[ndarray, ndarray] [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 low-pass 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 of shape (2M,) or of shape (2M,2M)) – covariance matrix or squared uncertainty associated with Y
UF (np.ndarray of shape (2M,2M), optional) – covariance matrix associated with F
- 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
- Raises:
ValueError – If dimensions of Y and F do not match.
- 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.GUM_DFT(x: ndarray, Ux: ndarray | float, N: int | None = None, window: ndarray | None = None, CxCos: ndarray | None = None, CxSin: ndarray | None = None, returnC: bool = False, mask: ndarray | None = None) Tuple[ndarray, Tuple[ndarray, ndarray, ndarray] | ndarray] | Tuple[ndarray, Tuple[ndarray, ndarray, ndarray] | ndarray, Dict[str, ndarray]] [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 (np.ndarray of shape (M,)) – vector of time domain signal values
Ux (np.ndarray of shape (M,) or of shape (M,M) or float) – covariance matrix associated with x, or vector of squared standard uncertainties, or noise variance as float
N (int, optional) – length of time domain signal for DFT; N>=len(x)
window (np.ndarray of shape (M,), optional) – vector of the time domain window values
CxCos (np.ndarray, optional) – cosine part of sensitivity matrix
CxSin (np.ndarray, optional) – sine part of sensitivity matrix
returnC (bool, optional) – if True, return sensitivity matrix blocks, if False (default) do not return them
mask (ndarray of dtype bool, optional) – calculate DFT values and uncertainties only at those frequencies where mask is True
- Returns:
F (np.ndarray) – vector of complex valued DFT values or of its real and imaginary parts
UF (np.ndarray) – covariance matrix associated with real and imaginary part of F
CxCos and CxSin (Dict) – Keys are “CxCos”, “CxSin” and values the respective sensitivity matrix entries
References
Eichstädt and Wilkens [Eichst2016]
- Raises:
ValueError – If N < len(x)
- PyDynamic.uncertainty.propagate_DFT.GUM_DFTfreq(N, dt=1)[source]
Return the Discrete Fourier Transform sample frequencies
- Parameters:
- Returns:
f – Array of length
n//2 + 1
containing the sample frequencies- Return type:
ndarray
See also
None
:numpy.fft.rfftfreq
- PyDynamic.uncertainty.propagate_DFT.GUM_iDFT(F: ndarray, UF: ndarray, Nx: int | None = None, Cc: ndarray | None = None, Cs: ndarray | None = None, returnC: bool = False) Tuple[ndarray, ndarray] | Tuple[ndarray, ndarray, Dict[str, ndarray]] [source]
Propagation of squared uncertainties UF associated with the DFT values F
GUM propagation of the squared uncertainties 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 real-valued 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) – length 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 (bool, optional) – If True, return sensitivity matrix blocks (without scaling factor 1/N), if False do not return them
- Returns:
x (np.ndarray) – vector of time domain signal values
Ux (np.ndarray) – covariance matrix associated with x
Cc and Cs (Dict) – Keys are “Cc”, “Cs” and values the respective sensitivity matrix entries
References
Eichstädt and Wilkens [Eichst2016]
- PyDynamic.uncertainty.propagate_DFT.Time2AmpPhase(x: ndarray, Ux: ndarray | float) Tuple[ndarray, ndarray, ndarray] [source]
Transformation from time domain to amplitude and phase via DFT
- Parameters:
x (np.ndarray of shape (N,)) – time domain signal
Ux (np.ndarray of shape (N,) or of shape (N,N) or float) – covariance matrix associated with x, or vector of squared standard uncertainties, or noise variance as float
- 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, nx)) – M time domain signals of length nx
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:N-1
- 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(UAA), diag(UAP), diag(UPP)]
Uncertainty evaluation for the DWT
This module assists in uncertainty propagation for the discrete wavelet transform
The PyDynamic.uncertainty.propagate_DWT
module implements methods for
the propagation of uncertainties in the application of the discrete wavelet
transform (DWT).
This modules contains the following functions:
dwt()
: single level DWTwave_dec()
: wavelet decomposition / multi level DWTwave_dec_realtime()
: multi level DWTinv_dwt()
: single level inverse DWTwave_rec()
: wavelet reconstruction / multi level inverse DWTfilter_design()
: provide common wavelet filters (viapywt.Wavelet
)dwt_max_level()
: return the maximum achievable DWT level
- PyDynamic.uncertainty.propagate_DWT.dwt(x, Ux, lowpass, highpass, states=None, realtime=False, subsample_start=1)[source]
Apply low-pass
lowpass
and high-passhighpass
to time-series datax
The uncertainty is propagated through the transformation by using
PyDynamic.uncertainty.propagate_filter.IIRuncFilter()
.Return the subsampled results.
- Parameters:
x (np.ndarray) – filter input signal
Ux (float or np.ndarray) – float: standard deviation of white noise in x 1D-array: point-wise standard uncertainties of non-stationary white noise
lowpass (np.ndarray) – FIR filter coefficients representing a low-pass for decomposition
highpass (np.ndarray) – FIR filter coefficients representing a high-pass for decomposition
states (dictionary of internal high/lowpass-filter states, optional (default: None)) – allows to continue at the last used internal state from previous call
realtime (Boolean, optional (default: False)) – for realtime applications, no signal padding has to be done before decomposition
subsample_start (int, optional (default: 1)) – At which position the subsampling should start, typically 1 (default) or 0. You should be happy with the default. We only need this to realize
wave_dec_realtime()
.
- Returns:
c_approx (np.ndarray) – subsampled low-pass output signal
U_approx (np.ndarray) – subsampled low-pass output uncertainty
c_detail (np.ndarray) – subsampled high-pass output signal
U_detail (np.ndarray) – subsampled high-pass output uncertainty
states (dictionary of internal high/lowpass-filter states) – allows to continue at the last used internal state in next call
- PyDynamic.uncertainty.propagate_DWT.dwt_max_level(data_length, filter_length)[source]
Return the highest achievable DWT level, given the provided data/filter lengths
- PyDynamic.uncertainty.propagate_DWT.filter_design(kind)[source]
Provide low- and highpass filters suitable for discrete wavelet transformation
This wraps
pywt.Wavelet
.- Parameters:
kind (string) –
filter name, i.e. db4, coif6, gaus9, rbio3.3, …
supported families:
pywt.families()
supported wavelets:
pywt.wavelist()
- Returns:
ld (np.ndarray) – low-pass filter for decomposition
hd (np.ndarray) – high-pass filter for decomposition
lr (np.ndarray) – low-pass filter for reconstruction
hr (np.ndarray) – high-pass filter for reconstruction
- PyDynamic.uncertainty.propagate_DWT.inv_dwt(c_approx, U_approx, c_detail, U_detail, lowpass, highpass, states=None, realtime=False)[source]
Single step of inverse discrete wavelet transform
- Parameters:
c_approx (np.ndarray) – low-pass output signal
U_approx (np.ndarray) – low-pass output uncertainty
c_detail (np.ndarray) – high-pass output signal
U_detail (np.ndarray) – high-pass output uncertainty
lowpass (np.ndarray) – FIR filter coefficients representing a low-pass for reconstruction
highpass (np.ndarray) – FIR filter coefficients representing a high-pass for reconstruction
states (dict, optional (default: None)) – internal high/lowpass-filter states, allows to continue at the last used internal state from previous call
realtime (Boolean, optional (default: False)) – for realtime applications, no signal padding has to be undone after reconstruction
- Returns:
x (np.ndarray) – upsampled reconstructed signal
Ux (np.ndarray) – upsampled uncertainty of reconstructed signal
states (dictionary of internal high/lowpass-filter states) – allows to continue at the last used internal state in next call
- PyDynamic.uncertainty.propagate_DWT.wave_dec(x, Ux, lowpass, highpass, n=-1)[source]
Multilevel discrete wavelet transformation of time-series x with uncertainty Ux
- Parameters:
x (np.ndarray) – input signal
Ux (float or np.ndarray) – float: standard deviation of white noise in x 1D-array: point-wise standard uncertainties of non-stationary white noise
lowpass (np.ndarray) – decomposition low-pass for wavelet_block
highpass (np.ndarray) – decomposition high-pass for wavelet_block
n (int, optional (default: -1)) – consecutive repetitions of wavelet_block user is warned, if it is not possible to reach the specified depth use highest possible level if set to -1 (default)
- Returns:
coeffs (list of arrays) – order of arrays within list is: [cAn, cDn, cDn-1, …, cD2, cD1]
Ucoeffs (list of arrays) – uncertainty of coeffs, same order as coeffs
original_length (int) – equals to len(x) necessary to restore correct length
- PyDynamic.uncertainty.propagate_DWT.wave_dec_realtime(x, Ux, lowpass, highpass, n=1, level_states=None)[source]
Multilevel discrete wavelet transformation of time-series x with uncertainty Ux
Similar to
wave_dec()
, but allows to start from the internal_state of a previous call.- Parameters:
x (np.ndarray) – input signal
Ux (float or np.ndarray) – float: standard deviation of white noise in x 1D-array: point-wise standard uncertainties of non-stationary white noise
lowpass (np.ndarray) – decomposition low-pass for wavelet_block
highpass (np.ndarray) – decomposition high-pass for wavelet_block
n (int, optional (default: 1)) – consecutive repetitions of wavelet_block There is no maximum level in continuous wavelet transform, so the default is n=1
level_states (dict, optional (default: None)) – internal state from previous call
- Returns:
coeffs (list of arrays) – order of arrays within list is: [cAn, cDn, cDn-1, …, cD2, cD1]
Ucoeffs (list of arrays) – uncertainty of coeffs, same order as coeffs
original_length (int) – equals to len(x) necessary to restore correct length
level_states (dict) – last internal state
- PyDynamic.uncertainty.propagate_DWT.wave_rec(coeffs, Ucoeffs, lowpass, highpass, original_length=None)[source]
Multilevel discrete wavelet reconstruction from levels back into time-series
- Parameters:
coeffs (list of arrays) –
order of arrays within list is: [cAn, cDn, cDn-1, …, cD2, cD1] where:
cAi: approximation coefficients array from i-th level
cDi: detail coefficients array from i-th level
Ucoeffs (list of arrays) – uncertainty of coeffs, same order as coeffs
lowpass (np.ndarray) – reconstruction low-pass for wavelet_block
highpass (np.ndarray) – reconstruction high-pass for wavelet_block
original_length (int, optional (default: None)) – necessary to restore correct length of original time-series
- Returns:
x (np.ndarray) – reconstructed signal
Ux (np.ndarray) – uncertainty of reconstructed signal
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 thetaIIRuncFilter()
: Uncertainty propagation for the signal x and the uncertain IIR filter (b,a)IIR_get_initial_state()
: Get a valid internal state forIIRuncFilter()
that assumes a stationary signal before the first value.
Note
The Elster-Link 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', return_full_covariance=False)[source]
Uncertainty propagation for signal y and uncertain FIR filter theta
A preceding FIR low-pass filter with coefficients blow can be provided optionally.
This method keeps the signature of PyDynamic.uncertainty.FIRuncFilter, but internally works differently and can return a full covariance matrix. Also, sigma_noise can be a full covariance matrix.
- Parameters:
y (np.ndarray) – filter input signal
sigma_noise (float or np.ndarray) –
float: standard deviation of white noise in y
1D-array: interpretation depends on kind
2D-array: full covariance of input
theta (np.ndarray) – FIR filter coefficients
Utheta (np.ndarray, optional) –
1D-array: coefficient-wise standard uncertainties of filter
2D-array: covariance matrix associated with theta
if the filter is fully certain, use Utheta = None (default) to make use of more efficient calculations. see also the comparison given in <examplesDigital filteringFIRuncFilter_runtime_comparison.py>
shift (int, optional) – time delay of filter output signal (in samples) (defaults to 0)
blow (np.ndarray, optional) – optional FIR low-pass filter
kind (string, optional) –
only meaningful in combination with sigma_noise a 1D numpy array
”diag”: point-wise standard uncertainties of non-stationary white noise
”corr”: single sided autocovariance of stationary colored noise (default)
return_full_covariance (bool, optional) – whether or not to return a full covariance of the output, defaults to False
- Returns:
x (np.ndarray) – FIR filter output signal
Ux (np.ndarray) –
return_full_covariance == False : point-wise standard uncertainties associated with x (default)
return_full_covariance == True : covariance matrix containing uncertainties associated with x
References
Elster and Link 2008 [Elster2008]
- PyDynamic.uncertainty.propagate_filter.IIR_get_initial_state(b, a, Uab=None, x0=1.0, U0=1.0, Ux=None)[source]
Calculate the internal state for the IIRuncFilter-function corresponding to stationary non-zero input signal.
- Parameters:
b (np.ndarray) – filter numerator coefficients
a (np.ndarray) – filter denominator coefficients
Uab (np.ndarray, optional (default: None)) – covariance matrix for (a[1:],b)
x0 (float, optional (default: 1.0)) – stationary input value
U0 (float, optional (default: 1.0)) – stationary input uncertainty
Ux (np.ndarray, optional (default: None)) – single sided autocovariance of stationary (colored/correlated) noise (needed in the kind=”corr” case of
IIRuncFilter()
)
- Returns:
internal_state – dictionary of state
- Return type:
- PyDynamic.uncertainty.propagate_filter.IIRuncFilter(x, Ux, b, a, Uab=None, state=None, kind='corr')[source]
Uncertainty propagation for the signal x and the uncertain IIR filter (b,a)
- Parameters:
x (np.ndarray) – filter input signal
Ux (float or np.ndarray) – float: standard deviation of white noise in x (requires kind=”diag”) 1D-array: interpretation depends on kind
b (np.ndarray) – filter numerator coefficients
a (np.ndarray) – filter denominator coefficients
Uab (np.ndarray, optional (default: None)) – covariance matrix for (a[1:],b)
state (dict, optional (default: None)) –
An internal state (z, dz, P, cache) to start from - e.g. from a previous run of IIRuncFilter.
If not given, (z, dz, P) are calculated such that the signal was constant before the given range
If given, the input parameters (b, a, Uab) are ignored to avoid repetitive rebuild of the internal system description (instead, the cache is used). However a valid new state (i.e. with new b, a, Uab) can always be generated by using
IIR_get_initial_state()
.
kind (string, optional (default: "corr")) – defines the interpretation of Ux, if Ux is a 1D-array “diag”: point-wise standard uncertainties of non-stationary white noise “corr”: single sided autocovariance of stationary (colored/correlated) noise (default)
- Returns:
y (np.ndarray) – filter output signal
Uy (np.ndarray) – uncertainty associated with y
state (dict) – dictionary of updated internal state
Note
In case of a == [1.0] (FIR filter), the results of
IIRuncFilter()
andFIRuncFilter()
might differ! This is because IIRuncFilter propagates uncertainty according to the (first-order Taylor series of the) GUM, whereas FIRuncFilter takes full variance information into account (which leads to an additional term). This is documented in the description of formula (33) of [Elster2008] . The difference can be visualized by running PyDynamic/examples/digital_filtering/validate_FIR_IIR_MC.pyReferences
Link and Elster [Link2009]
Monte Carlo methods for digital filtering
“Monte Carlo methods for the propagation of uncertainties 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 GUM-S2 Monte Carlo method provides a reference method for the calculation of uncertainties for such cases.
This module contains the following functions:
MC()
: Standard Monte Carlo method for application of digital filterSMC()
: Sequential Monte Carlo method with reduced computer memory requirementsUMC()
: Update Monte Carlo method for application of digital filters with reduced computer memory requirementsUMC_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), point-wise 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
- Returns:
y, Uy (np.ndarray) – filtered output signal and associated uncertainties, only returned if return_samples is
False
Y (np.ndarray) – array of Monte Carlo results, only returned if return_samples is
True
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 low-pass filter numerator coefficients
alow (np.ndarray) – optional low-pass 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
phi (np.ndarray, optional) – parameters for AR(MA) noise model \(\epsilon(n) = \sum_k \phi_k\epsilon(n-k) + \sum_k \theta_k w(n-k) + w(n)\) with \(w(n)\sim N(0,noise_std^2)\)
theta (np.ndarray, optional) – parameters for AR(MA) noise model \(\epsilon(n) = \sum_k \phi_k\epsilon(n-k) + \sum_k \theta_k w(n-k) + 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 point-wise 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: bin-edges val[“bin-edges”], bin-counts val[“bin-counts”]
References
Eichstädt, Link, Harris, Elster [Eichst2012]
ported to python in 2019-08 from matlab-version of Sascha Eichstaedt (PTB) from 2011-10-12
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=2, return_histograms=True, compute_full_covariance=True)[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 multi-dimensional, 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 return-value of documentation
n_cpu (int, optional) – number of CPUs to use for multiprocessing, defaults to all available CPUs
return_histograms (bool, optional) – whether to compute a histogram for each entry of the result at all
compute_full_covariance (bool, optional) – whether to compute the full covariance matrix or just its diagonal
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 bin-edges and bin-counts
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 memory-improvements of the UMC-method.- Returns:
sims – dict of samples and corresponding results of every evaluated simulation samples and results are saved in their original shape
- Return type:
References
Eichstädt, Link, Harris, Elster [Eichst2012]
Uncertainty evaluation for interpolation
This module assists in uncertainty propagation for 1-dimensional interpolation
The PyDynamic.uncertainty.interpolate
module implements methods for
the propagation of uncertainties in the application of standard interpolation methods
as provided by scipy.interpolate.interp1d
.
This module contains the following functions:
interp1d_unc()
: Interpolate arbitrary time series considering the associated uncertaintiesmake_equidistant()
: Interpolate a 1-D function equidistantly considering associated uncertainties
- PyDynamic.uncertainty.interpolate.interp1d_unc(x_new: ndarray, x: ndarray, y: ndarray, uy: ndarray, kind: str | None = 'linear', copy=True, bounds_error: bool | None = None, fill_value: float | Tuple[float, float] | str | None = nan, fill_unc: float | Tuple[float, float] | str | None = nan, assume_sorted: bool | None = True, returnC: bool | None = False) Tuple[ndarray, ndarray, ndarray] | Tuple[ndarray, ndarray, ndarray, ndarray] [source]
Interpolate a 1-D function considering the associated uncertainties
x and y are arrays of values used to approximate some function \(f \colon y = f(x)\).
Note that calling
interp1d_unc()
with NaNs present in input values results in undefined behaviour.An equal number of each of the original x and y values and associated uncertainties is required.
- Parameters:
x_new ((M,) array_like) – A 1-D array of real values to evaluate the interpolant at. x_new can be sorted in any order.
x ((N,) array_like) – A 1-D array of real values.
y ((N,) array_like) – A 1-D array of real values. The length of y must be equal to the length of x.
uy ((N,) array_like) – A 1-D array of real values representing the standard uncertainties associated with y.
kind (str, optional) – Specifies the kind of interpolation for y as a string (‘previous’, ‘next’, ‘nearest’, ‘linear’ or ‘cubic’). Default is ‘linear’.
copy (bool, optional) – If True, the method makes internal copies of x and y. If False, references to x and y are used. The default is to copy.
bounds_error (bool, optional) – If True, a ValueError is raised any time interpolation is attempted on a value outside of the range of x (where extrapolation is necessary). If False, out of bounds values are assigned fill_value. By default, an error is raised unless
fill_value="extrapolate"
.fill_value (array-like or (array-like, array_like) or “extrapolate”, optional) –
if a ndarray (or float), this value will be used to fill in for requested points outside of the data range. If not provided, then the default is NaN. The array-like must broadcast properly to the dimensions of the non-interpolation axes.
If a two-element tuple, then the first element is used as a fill value for
x_new < t[0]
and the second element is used forx_new > t[-1]
. Anything that is not a 2-element tuple (e.g., list or ndarray, regardless of shape) is taken to be a single array-like argument meant to be used for both bounds asbelow, above = fill_value, fill_value
.If “extrapolate”, then points outside the data range will be set to the first or last element of the values.
If cubic-interpolation, C2-continuity at the transition to the extrapolation-range is not guaranteed. This behavior might change in future implementations, see issue #210 for details.
Both parameters fill_value and fill_unc should be provided to ensure desired behaviour in the extrapolation range.
fill_unc (array-like or (array-like, array_like) or “extrapolate”, optional) – Usage and behaviour as described in fill_value but for the uncertainties. Both parameters fill_value and fill_unc should be provided to ensure desired behaviour in the extrapolation range.
assume_sorted (bool, optional) – If False, values of x can be in any order and they are sorted first. If True, x has to be an array of monotonically increasing values.
returnC (bool, optional) – If True, return sensitivity coefficients for later use. This is only available for interpolation kind ‘linear’ and for fill_unc=”extrapolate” at the moment. If False sensitivity coefficients are not returned and internal computation is slightly more efficient.
- Returns:
x_new ((M,) array_like) – values at which the interpolant is evaluated
y_new ((M,) array_like) – interpolated values
uy_new ((M,) array_like) – interpolated associated standard uncertainties
C ((M,N) array_like) – sensitivity matrix \(C\), which is used to compute the uncertainties \(U_{y_{new}} = C \cdot \operatorname{diag}(u_y^2) \cdot C^T\), only returned if returnC is False, which is the default behaviour.
References
White [White2017]
- PyDynamic.uncertainty.interpolate.make_equidistant(x: ndarray, y: ndarray, uy: ndarray, dx: float | None = 0.05, kind: str | None = 'linear') Tuple[ndarray, ndarray, ndarray] [source]
Interpolate a 1-D function equidistantly considering associated uncertainties
Interpolate function values equidistantly and propagate uncertainties accordingly.
x and y are arrays of values used to approximate some function \(f \colon y = f(x)\).
Note that calling
interp1d_unc()
with NaNs present in input values results in undefined behaviour.An equal number of each of the original x and y values and associated uncertainties is required.
- Parameters:
x ((N,) array_like) – A 1-D array of real values.
y ((N,) array_like) – A 1-D array of real values. The length of y must be equal to the length of x.
uy ((N,) array_like) – A 1-D array of real values representing the standard uncertainties associated with y.
dx (float, optional) – desired interval length (defaults to 5e-2)
kind (str, optional) – Specifies the kind of interpolation for y as a string (‘previous’, ‘next’, ‘nearest’, ‘linear’ or ‘cubic’). Default is ‘linear’.
- Returns:
x_new ((M,) array_like) – values at which the interpolant is evaluated
y_new ((M,) array_like) – interpolated values
uy_new ((M,) array_like) – interpolated associated standard uncertainties
References
White [White2017]