# -*- coding: utf-8 -*-
"""
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 implements Monte Carlo methods for the propagation of uncertainties for digital filtering.
"""
# TODO: Implement updating Monte Carlo method
<<<<<<< HEAD
<<<<<<< Updated upstream
# TODO: Implement repeated random number drawing from multivariate normal like in mvnrnd with re-using cholesky
=======
# TODO: Advocate use of second-order-systems for higher-order IIR filters (i.e. throw warning at the user)
>>>>>>> devel1
import numpy as np
import scipy as sp
=======
import numpy as np
import scipy as sp
from numpy import matrix
>>>>>>> Stashed changes
import sys
from scipy.signal import lfilter
import scipy.stats as stats
from ..misc.tools import zerom
from ..misc.filterstuff import isstable
__all__ = ["MC", "SMC"]
class Normal_ZeroCorr():
"""
Multivariate normal distribution with zero correlation
"""
<<<<<<< Updated upstream
def __init__(self, mean=None, cov=None):
=======
def __init__(self, loc=None, scale=None):
>>>>>>> Stashed changes
"""
Parameters
----------
loc: np.ndarray, optional
mean values, default is zero
scale: np.ndarray, optional
standard deviations for the elements in loc, default is zero
"""
if isinstance(mean, np.ndarray):
self.mean = mean
if isinstance(cov, np.ndarray):
assert (len(cov)==len(mean))
self.std = cov
elif isinstance(cov, float):
self.std = cov * np.ones_like(mean)
else:
self.std = np.zeros_like(mean)
elif isinstance(cov, np.ndarray):
self.std = cov
self.mean = np.zeros_like(cov)
def rvs(self, size=1):
return np.tile(self.mean, (size, 1)) + np.random.randn(size, len(self.mean))*np.tile(self.std, (size, 1))
class Normal_FullCov(stats._multivariate.multivariate_normal_frozen):
"""
Derive from the default multivariate normal class one that allows storing the covariance matrix decomposition once calculated
"""
def __init__(self, mean=0, cov=1, seed=None, allow_singular=False):
super(Normal_FullCov, self).__init__(mean, cov, allow_singular, seed)
# initialize covariance matrix decomposition
[docs]def MC(x,Ux,b,a,Uab,runs=1000,blow=None,alow=None,return_samples=False,shift=0,verbose=True):
"""Standard Monte Carlo method
Monte Carlo based propagation of uncertainties for a digital filter (b,a)
with uncertainty matrix
:math:`U_{\theta}` for :math:`\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 :math:`U_\theta`
runs: int,optional
number of Monte Carlo runs
return_samples: bool, optional
whether samples or mean and std are returned
If 'return_sampes' is false, the method returns:
Returns
-------
y: np.ndarray
filter output signal
Uy: np.ndarray
uncertainty associated with
Other wise the method returns
Returns
-------
Y: np.ndarray
array of Monte Carlo results
References
----------
* Eichstädt, Link, Harris and Elster [Eichst2012]_
"""
Na = len(a)
runs = int(runs)
Y = np.zeros((runs,len(x)))
theta = np.hstack((a[1:],b))
Theta = np.random.multivariate_normal(theta,Uab,runs)
if isinstance(Ux, np.ndarray):
<<<<<<< HEAD
<<<<<<< Updated upstream
dist = stats.multivariate_normal(x, Ux, allow_singular=True)
=======
=======
>>>>>>> devel1
if len(Ux.shape)==1:
dist = Normal_ZeroCorr(loc=x, scale=Ux)
else:
dist = stats.multivariate_normal(x, Ux)
<<<<<<< HEAD
>>>>>>> Stashed changes
=======
>>>>>>> devel1
elif isinstance(Ux, float):
dist = Normal_ZeroCorr(loc=x, scale=Ux)
else:
raise NotImplementedError("The supplied type of uncertainty is not implemented")
unst_count = 0
st_inds = list()
if verbose:
sys.stdout.write('MC progress: ')
for k in range(runs):
xn = dist.rvs()
if not blow is None:
if alow is None:
alow = 1.0
xn = lfilter(blow,alow,xn)
bb = Theta[k,Na-1:]
aa = np.hstack((1.0, Theta[k,:Na-1]))
if isstable(bb,aa):
Y[k,:] = lfilter(bb,aa,xn)
st_inds.append(k)
else:
unst_count += 1
if np.mod(k, 0.1*runs) == 0 and verbose:
sys.stdout.write(' %d%%' % (np.round(100.0*k/runs)))
if verbose:
sys.stdout.write(" 100%\n")
if unst_count > 0:
print("In %d Monte Carlo %d filters have been unstable" % (runs,unst_count))
print("These results will not be considered for calculation of mean and std")
print("However, if return_samples is 'True' then ALL samples are returned.")
Y = np.roll(Y,int(shift),axis=1)
if return_samples:
return Y
else:
y = np.mean(Y[st_inds,:],axis=0)
uy= np.std(Y[st_inds,:],axis=0)
return y, uy
[docs]def 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):
"""Sequential Monte Carlo method
Sequential Monte Carlo propagation for a digital filter (b,a) with uncertainty
matrix :math:`U_{\theta}` for :math:`\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 :math:`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, theta: np.ndarray, optional
parameters for AR(MA) noise model
::math:`\epsilon(n) = \sum_k \phi_k\epsilon(n-k) + \sum_k \theta_k w(n-k) + w(n)`
with ::math:`w(n)\sim N(0,noise_std^2)`
Delta: float,optional
upper bound on systematic error of the filter
If return_samples is False:
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 not None) at
Otherwise:
Returns
-------
Y: np.ndarray
array of all Monte Carlo results
References
----------
* Eichstädt, Link, Harris, Elster [Eichst2012]_
"""
runs = int(runs)
if isinstance(a,np.ndarray):
Na = len(a)-1
else:
Na = 0
if isinstance(b,np.ndarray):
Nb = len(b)-1
else:
Nb = 0
if isinstance(theta,np.ndarray) or isinstance(theta,float):
if isinstance(theta,float):
W = np.zeros((runs,1))
else:
W = np.zeros((runs,len(theta)))
else:
MA = False
if isinstance(phi,np.ndarray) or isinstance(phi,float):
AR = True
if isinstance(phi,float):
E = np.zeros((runs,1))
else:
E = np.zeros((runs,len(phi)))
else:
AR = False
if isinstance(blow,np.ndarray):
X = np.zeros((runs,len(blow)))
else:
X = np.zeros(runs)
if isinstance(alow,np.ndarray):
Xl = np.zeros((runs,len(alow)-1))
else:
Xl = np.zeros((runs,1))
if Na==0:
coefs = b
else:
coefs = np.hstack((a[1:],b))
if isinstance(Uab, np.ndarray):
Coefs = np.random.multivariate_normal(coefs, Uab, runs)
else:
Coefs = np.tile(coefs, (runs, 1))
b0 = Coefs[:,Na]
if Na>0: # filter is IIR
A = Coefs[:,:Na]
if Nb>Na:
A = np.hstack((A,zerom((runs,Nb-Na))))
else: # filter is FIR -> zero state equations
A = np.zeros((runs,Nb))
c = Coefs[:,Na+1:] - np.multiply(np.tile(b0[:, np.newaxis],(1,Nb)),A)
States = np.zeros(np.shape(A))
calcP = False
if not Perc is None:
calcP = True
P = np.zeros((len(Perc),len(x)))
y = np.zeros_like(x)
Uy= np.zeros_like(x)
print("Sequential Monte Carlo progress", end="")
for k in range(len(x)):
w = np.random.randn(runs)*noise_std
if AR and MA:
E = np.hstack( ( E.dot(phi) + W.dot(theta) + w, E[:-1]) )
W = np.hstack( (w, W[:-1] ) )
elif AR:
E = np.hstack( (E.dot(phi) + w, E[:-1]) )
elif MA:
E = W.dot(theta) + w
W = np.hstack( (w, W[:-1] ) )
else:
w = np.random.randn(runs,1)*noise_std
E = w
if isinstance(alow,np.ndarray):
X = np.hstack((x[k] + E, X[:, :-1]))
Xl = np.hstack( ( X.dot(blow.T) - Xl[:,:len(alow)].dot(alow[1:]), Xl[:,:-1] ) )
elif isinstance(blow,np.ndarray):
X = np.hstack((x[k] + E, X[:, :-1]))
Xl = X.dot(blow.T)
else:
Xl = x[k] + E
Y = np.sum(np.multiply(c,States),axis=1) + np.multiply(b0,Xl[:,0]) + (np.random.rand(runs)*2*Delta - Delta)
Z = -np.sum(np.multiply(A,States),axis=1) + Xl[:,0]
States = np.hstack((Z[:,np.newaxis], States[:,:-1]))
y[k] = np.mean(Y)
Uy[k]= np.std(Y)
if calcP:
P[:,k] = sp.stats.mstats.mquantiles(np.asarray(Y),prob=Perc)
if np.mod(k, np.round(0.1*len(x))) == 0:
print(' %d%%' % (np.round(100.0*k/len(x))), end="")
print(" 100%")
y = np.roll(y,int(shift))
Uy= np.roll(Uy,int(shift))
if calcP:
P = np.roll(P,int(shift),axis=1)
return y, Uy, P
else:
return y, Uy