Miscellaneous utilities for time series analysis.
Perform an iterative procedure to find the optimal weights for K direct spectral estimators of DPSS tapered signals.
Parameters : | yk : ndarray (K, N)
eigvals : ndarray, length-K
sides : str
max_iter : int
|
---|---|
Returns : | weights, nu :
|
Notes
The weights to use for making the multitaper estimate, such that S_{mt} = \sum_{k} |w_k|^2S_k^{mt} / \sum_{k} |w_k|^2
If there are less than 3 tapers, then the adaptive weights are not found. The square root of the eigenvalues are returned as weights, and the degrees of freedom are 2*K
A measure of the goodness of fit of an auto-regressive model based on the model order and the error covariance.
Parameters : | ecov : float array
p :
m : int
Ntotal :
corrected : boolean (optional)
|
---|---|
Returns : | AIC : float
|
Notes
This is an implementation of equation (50) in Ding et al. (2006):
M Ding and Y Chen and S Bressler (2006) Granger Causality: Basic Theory and Application to Neuroscience. http://arxiv.org/abs/q-bio/0608035v1
Correction for small sample size is taken from: http://en.wikipedia.org/wiki/Akaike_information_criterion.
Make an anti-symmetric random 2-d array of shape (size,size).
Parameters : | n : int
sample_func : function, optional.
|
---|
Examples
>>> np.random.seed(0) # for doctesting
>>> np.set_printoptions(precision=4) # for doctesting
>>> antisymm_rand_arr(4)
array([[ 0. , 0.7152, 0.6028, 0.5449],
[-0.7152, 0. , 0.4376, 0.8918],
[-0.6028, -0.4376, 0. , 0.5289],
[-0.5449, -0.8918, -0.5289, 0. ]])
This generates a signal u(n) = a1*u(n-1) + a2*u(n-2) + ... + v(n) where v(n) is a stationary stochastic process with zero mean and variance = sigma. XXX: confusing variance notation
Parameters : | N : float
sigma : float
coefs : list or array of floats
drop_transients : float
v : float array
|
---|---|
Returns : | u : ndarray
v : ndarray
coefs : ndarray
The form of the feedback coefficients is a little different than : the normal linear constant-coefficient difference equation. Therefore : the transfer function implemented in this method is : H(z) = sigma**0.5 / ( 1 - sum_k coefs(k)z**(-k) ) 1 <= k <= P : |
Examples
>>> import nitime.algorithms as alg
>>> ar_seq, nz, alpha = ar_generator()
>>> fgrid, hz = alg.freq_response(1.0, a=np.r_[1, -alpha])
>>> sdf_ar = (hz * hz.conj()).real
Returns the autocorrelation of signal s at all lags.
Parameters : | x : ndarray axis : time axis all_lags : {True/False}
|
---|
Notes
Adheres to the definition
R_{xx}[k]=E{X[n+k]X^{*}[n]}
where X is a discrete, stationary (ergodic) random process
Returns the autocovariance of signal s at all lags.
Parameters : | x : ndarray axis : time axis all_lags : {True/False}
|
---|---|
Returns : | cxx : ndarray
|
Notes
Adheres to the definition
C_{xx}[k]=E{(X[n+k]-E{X})(X[n]-E{X})^{*}}
where X is a discrete, stationary (ergodic) random process
This method computes the following function
R_{xx}(k) = E{ x(t)x^{*}(t-k) } = E{ x(t+k)x^{*}(t) } k in {0, 1, ..., nlags-1}
(* := conjugate transpose)
Note: this is related to the other commonly used definition for vector autocovariance
R_{xx}^{(2)}(k) = E{ x(t-k)x^{*}(t) } = R_{xx}(-k) = R_{xx}^{*}(k)
Parameters : | x : ndarray (nc, N) nlags : int, optional
|
---|---|
Returns : | rxx : ndarray (nc, nc, nlags) |
The Bayesian Information Criterion, also known as the Schwarz criterion is a measure of goodness of fit of a statistical model, based on the
number of model parameters and the likelihood of the model
Parameters : | ecov : float array
|
---|---|
Returns : | BIC : float
|
rac{2p^2 m log(N_{total})}{N_{total}},
where $Sigma$ is the noise covariance matrix. In auto-regressive model estimation, this matrix will contain in $Sigma_{i,j}$ the residual variance in estimating time-series $i$ from $j$, $p$ is the dimensionality of the data, $m$ is the number of parameters in the model and $N_{total}$ is the number of time-points.
M Ding and Y Chen and S Bressler (2006) Granger Causality: Basic Theory and Application to Neuroscience. http://arxiv.org/abs/q-bio/0608035v1
For a frequency grid spaced on the unit circle of an imaginary plane, return the corresponding freqency grid in Hz.
Maps the input into the continuous interval (bottom, top) where bottom defaults to 0 and top defaults to 2*pi
Parameters : | x : ndarray - the input array bottom : float, optional (defaults to 0).
top : float, optional (defaults to 2*pi).
|
---|---|
Returns : | The input array, mapped into the interval (bottom,top) : |
Returns the crosscorrelation sequence between two ndarrays. This is performed by calling fftconvolve on x, y[::-1]
Parameters : | x : ndarray y : ndarray axis : time axis all_lags : {True/False}
|
---|---|
Returns : | rxy : ndarray
|
Notes
cross correlation is defined as
R_{xy}[k]=E{X[n+k]Y^{*}[n]}
where X and Y are discrete, stationary (ergodic) random processes
Returns the crosscovariance sequence between two ndarrays. This is performed by calling fftconvolve on x, y[::-1]
Parameters : | x : ndarray y : ndarray axis : time axis all_lags : {True/False}
debias : {True/False}
|
---|---|
Returns : | cxy : ndarray
|
Notes
cross covariance of processes x and y is defined as
C_{xy}[k]=E{(X(n+k)-E{X})(Y(n)-E{Y})^{*}}
where X and Y are discrete, stationary (or ergodic) random processes
Also note that this routine is the workhorse for all auto/cross/cov/corr functions.
This method computes the following function
R_{xy}(k) = E{ x(t)y^{*}(t-k) } = E{ x(t+k)y^{*}(t) } k \in {0, 1, ..., nlags-1}
(* := conjugate transpose)
Note: This is related to the other commonly used definition for vector crosscovariance
R_{xy}^{(2)}(k) = E{ x(t-k)y^{*}(t) } = R_{xy}^(-k) = R_{yx}^{*}(k)
Parameters : | x, y : ndarray (nc, N) nlags : int, optional
|
---|---|
Returns : | rxy : ndarray (nc, nc, nlags) |
Convert the values in x to decibels. If the values in x are in ‘power’-like units, then set the power flag accordingly
Detect the presence of line spectra in s using the F-test described in “Spectrum estimation and harmonic analysis” (Thompson 81). Strategies for detecting harmonics in low SNR include increasing the number of FFT points (NFFT keyword arg) and/or increasing the stability of the spectral estimate by using more tapers (higher NW parameter).
Returns : | (freq, beta) : sequence
|
---|
Return the indices to access the main diagonal of an array.
This returns a tuple of indices that can be used to access the main diagonal of an array with ndim (>=2) dimensions and shape (n,n,...,n). For ndim=2 this is the usual diagonal, for ndim>2 this is the set of indices to access A[i,i,...,i] for i=[0..n-1].
Parameters : | n : int
ndim : int, optional
|
---|
See also
-, array.
Examples
Create a set of indices to access the diagonal of a (4,4) array: >>> di = diag_indices(4)
>>> a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
>>> a
array([[ 1, 2, 3, 4],
[ 5, 6, 7, 8],
[ 9, 10, 11, 12],
[13, 14, 15, 16]])
>>> a[di] = 100
>>> a
array([[100, 2, 3, 4],
[ 5, 100, 7, 8],
[ 9, 10, 100, 12],
[ 13, 14, 15, 100]])
Now, we create indices to manipulate a 3-d array: >>> d3 = diag_indices(2,3)
And use it to set the diagonal of a zeros array to 1: >>> a = np.zeros((2,2,2),int) >>> a[d3] = 1 >>> a array([[[1, 0],
[0, 0]],
Return the indices to access the main diagonal of an n-dimensional array.
See diag_indices() for full details.
Parameters : | arr : array, at least 2-d |
---|
Compute the expected value of the jackknife variance estimate over K windows below. This expected value formula is based on the asymptotic expansion of the trigamma function derived in [Thompson_1994]
Returns : | evar : float
|
---|
Convolve two N-dimensional arrays using FFT. See convolve.
This is a fix of scipy.signal.fftconvolve, adding an axis argument and importing locally the stuff only needed for this function
Fill the main diagonal of the given array of any dimensionality.
For an array with ndim > 2, the diagonal is the list of locations with indices a[i,i,...,i], all identical.
This function modifies the input array in-place, it does not return a value.
This functionality can be obtained via diag_indices(), but internally this version uses a much faster implementation that never constructs the indices and uses simple slicing.
Parameters : | a : array, at least 2-dimensional.
val : scalar
|
---|
See also
-, -
Examples
>>> a = np.zeros((3,3),int)
>>> fill_diagonal(a,5)
>>> a
array([[5, 0, 0],
[0, 5, 0],
[0, 0, 5]])
The same function can operate on a 4-d array: >>> a = np.zeros((3,3,3,3),int) >>> fill_diagonal(a,4)
We only show a few blocks for clarity: >>> a[0,0] array([[4, 0, 0],
[0, 0, 0], [0, 0, 0]])
>>> a[1,1]
array([[0, 0, 0],
[0, 4, 0],
[0, 0, 0]])
>>> a[2,2]
array([[0, 0, 0],
[0, 0, 0],
[0, 0, 4]])
Create a FIR event matrix from a time-series of events.
Parameters : | events : 1-d int array
len_hrf : int
|
---|---|
Returns : | fir_matrix : matrix
|
Generates a multivariate autoregressive dataset given the formula:
X(t) + sum_{i=1}^{P} a(i)X(t-i) = E(t)
Where E(t) is a vector of samples from possibly covarying noise processes.
Parameters : | a : ndarray (n_order, n_c, n_c)
cov : ndarray (n_c, n_c)
N : int
|
---|---|
Returns : | mar, nz : mar and noise process shaped (n_c, N) : |
Find the indices of the lower and upper bounds within an array f
Parameters : | f, array : lb,ub, float : |
---|---|
Returns : | lb_idx, ub_idx: the indices into ‘f’ which correspond to values bounded : between ub and lb in that array : |
Returns the center frequencies of the frequency decomposotion of a time series of length n, sampled at Fs Hz
Calculate the analytical spectrum of a Hanning window
Parameters : | N : int
Fs : float
|
---|---|
Returns : | float array - the frequency bands, given N and FS : complex array: the power in the spectrum of the square window in the : frequency bands : |
Notes
This is equation 28b in Harris (1978):
W(\theta) = 0.5 D(\theta) + 0.25 (D(\theta - \frac{2\pi}{N}) + D(\theta + \frac{2\pi}{N}) ),
where:
D(\theta) = exp(j\frac{\theta}{2}) \frac{sin\frac{N\theta}{2}}{sin\frac{\theta}{2}}
F.J. Harris (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, 66:51-83
This is a verbatim copy of scipy.signal.hilbert from scipy version 0.8dev, which we carry around in order to use in case the version of scipy installed is old enough to have a broken implementation of hilbert
For two sets of coordinates, find the coordinates that are common to both, where the dimensionality is the coords1.shape[0]
Returns the variance of the coherency between x and y, estimated through jack-knifing the tapered samples in {tx, ty}.
Parameters : | tx : ndarray, (K, L)
ty : ndarray, (K, L)
eigvals : ndarray (K,)
|
---|---|
Returns : | jk_var : ndarray
|
Returns the variance of the log-sdf estimated through jack-knifing a group of independent sdf estimates.
Parameters : | yk : ndarray (K, L)
eigvals : ndarray (K,)
sides : str, optional
adpative : bool, optional
|
---|---|
Returns : | var : The estimate for log-sdf variance |
Notes
The jackknifed mean estimate is distributed about the true mean as a Student’s t-distribution with (K-1) degrees of freedom, and standard error equal to sqrt(var). However, Thompson and Chave [1] point out that this variance better describes the sample mean.
[1] Thomson D J, Chave A D (1991) Advances in Spectrum Analysis and Array Processing (Prentice-Hall, Englewood Cliffs, NJ), 1, pp 58-113.
Return the indices to access (n,n) arrays, given a masking function.
Assume mask_func() is a function that, for a square array a of size (n,n) with a possible offset argument k, when called as mask_func(a,k) returns a new array with zeros in certain locations (functions like triu() or tril() do precisely this). Then this function returns the indices where the non-zero values would be located.
Parameters : | n : int
mask_func : callable
k : scalar
|
---|---|
Returns : | indices : an n-tuple of index arrays.
|
Examples
These are the indices that would allow you to access the upper triangular part of any 3x3 array: >>> iu = mask_indices(3,np.triu)
For example, if a is a 3x3 array: >>> a = np.arange(9).reshape(3,3) >>> a array([[0, 1, 2],
[3, 4, 5], [6, 7, 8]])
Then: >>> a[iu] array([0, 1, 2, 4, 5, 8])
An offset can be passed also to the masking function. This gets us the indices starting on the first diagonal right of the main one: >>> iu1 = mask_indices(3,np.triu,1)
with which we now extract only three elements: >>> a[iu1] array([1, 2, 5])
Minmax_norm an array to [0,1] range.
By default, this simply rescales the input array to [0,1]. But it has a special ‘folding’ mode that allows for the normalization of an array with negative and positive values by mapping the negative values to their flipped sign
Parameters : | arr : 1d array mode : string, one of [‘direct’,’folding’] folding_edges : (float,float)
|
---|
Examples
>>> np.set_printoptions(precision=4) # for doctesting
>>> a = np.linspace(0.3,0.8,4)
>>> minmax_norm(a)
array([ 0. , 0.3333, 0.6667, 1. ])
>>> b = np.concatenate([np.linspace(-0.7,-0.3,3),
... np.linspace(0.3,0.8,3)])
>>> b
array([-0.7 , -0.5 , -0.3 , 0.3 , 0.55, 0.8 ])
>>> minmax_norm(b,'folding',[-0.3,0.3])
array([ 0.8, 0.4, 0. , 0. , 0.5, 1. ])
A function for finding the intersection of several different arrays
Parameters : | input is a tuple of arrays, with all the different arrays : |
---|---|
Returns : | array - the intersection of the inputs : |
Notes
Simply runs intersect1d iteratively on the inputs
The inverse transform of the above normalization
The generally accepted choice to transform coherence measures into a more normal distribution
Parameters : | x : ndarray, real
dof : int
copy : bool
|
---|---|
Returns : | y : ndarray, real
|
Returns the % signal change of each point of the times series along a given axis of the array time_series
Parameters : | ts : ndarray
ax : int, optional (default to -1)
|
---|---|
Returns : | ndarray :
|
Examples
>>> ts = np.arange(4*5).reshape(4,5)
>>> ax = 0
>>> percent_change(ts,ax)
array([[-100. , -88.2353, -78.9474, -71.4286, -65.2174],
[ -33.3333, -29.4118, -26.3158, -23.8095, -21.7391],
[ 33.3333, 29.4118, 26.3158, 23.8095, 21.7391],
[ 100. , 88.2353, 78.9474, 71.4286, 65.2174]])
>>> ax = 1
>>> percent_change(ts,ax)
array([[-100. , -50. , 0. , 50. , 100. ],
[ -28.5714, -14.2857, 0. , 14.2857, 28.5714],
[ -16.6667, -8.3333, 0. , 8.3333, 16.6667],
[ -11.7647, -5.8824, 0. , 5.8824, 11.7647]])
Subtracts an estimate of the mean from signal x at axis
Rescale an array to a new range.
Return a new array whose range of values is (amin,amax).
Parameters : | arr : array-like amin : float
amax : float
|
---|
Examples
>>> a = np.arange(5)
>>> rescale_arr(a,3,6)
array([ 3. , 3.75, 4.5 , 5.25, 6. ])
Calculate the analytical spectrum of a square window
Parameters : | N : int
Fs : float
|
---|---|
Returns : | float array - the frequency bands, given N and FS : complex array: the power in the spectrum of the square window in the : frequency bands : |
Notes
This is equation 21c in Harris (1978):
W(\theta) = exp(-j \frac{N-1}{2} \theta) \frac{sin \frac{N\theta}{2}} {sin\frac{\theta}{2}}
F.J. Harris (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, 66:51-83
Make a structured random 2-d array of shape (size,size).
If no optional arguments are given, a symmetric array is returned.
Parameters : | size : int
sample_func : function, optional.
utfac : float, optional
ltfac : float, optional
fill_diag : float, optional
|
---|
Examples
>>> np.random.seed(0) # for doctesting
>>> np.set_printoptions(precision=4) # for doctesting
>>> structured_rand_arr(4)
array([[ 0.5488, 0.7152, 0.6028, 0.5449],
[ 0.7152, 0.6459, 0.4376, 0.8918],
[ 0.6028, 0.4376, 0.7917, 0.5289],
[ 0.5449, 0.8918, 0.5289, 0.0871]])
>>> structured_rand_arr(4,ltfac=-10,utfac=10,fill_diag=0.5)
array([[ 0.5 , 8.3262, 7.7816, 8.7001],
[-8.3262, 0.5 , 4.6148, 7.8053],
[-7.7816, -4.6148, 0.5 , 9.4467],
[-8.7001, -7.8053, -9.4467, 0.5 ]])
Make a symmetric random 2-d array of shape (size,size).
Parameters : | n : int
sample_func : function, optional.
fill_diag : float, optional
|
---|
Examples
>>> np.random.seed(0) # for doctesting
>>> np.set_printoptions(precision=4) # for doctesting
>>> symm_rand_arr(4)
array([[ 0.5488, 0.7152, 0.6028, 0.5449],
[ 0.7152, 0.6459, 0.4376, 0.8918],
[ 0.6028, 0.4376, 0.7917, 0.5289],
[ 0.5449, 0.8918, 0.5289, 0.0871]])
>>> symm_rand_arr(4,fill_diag=4)
array([[ 4. , 0.8326, 0.7782, 0.87 ],
[ 0.8326, 4. , 0.4615, 0.7805],
[ 0.7782, 0.4615, 4. , 0.9447],
[ 0.87 , 0.7805, 0.9447, 4. ]])
Threshold values from the input array.
Parameters : | cmat : array threshold : float, optional.
threshold2 : float, optional.
|
---|---|
Returns : | indices, values: a tuple with ndim+1 : |
Examples
>>> np.set_printoptions(precision=4) # For doctesting
>>> a = np.linspace(0,0.2,5)
>>> a
array([ 0. , 0.05, 0.1 , 0.15, 0.2 ])
>>> threshold_arr(a,0.1)
(array([3, 4]), array([ 0.15, 0.2 ]))
With two thresholds: >>> threshold_arr(a,0.1,0.2) (array([0, 1]), array([ 0. , 0.05]))
Threshold values from the input matrix and return a new matrix.
Parameters : | arr : array threshold : float
threshold2 : float, optional.
|
---|---|
Returns : | An array shaped like the input, with the values outside the threshold : replaced with fill_val. : |
Perform an inverse iteration to find the eigenvector corresponding to the given eigenvalue in a symmetric tridiagonal system.
Parameters : | d : ndarray
e : ndarray
w : float
x0 : ndarray
rtol : float
|
---|---|
Returns : | e : ndarray
|
Return the indices for the lower-triangle of an (n,n) array.
Parameters : | n : int
k : int, optional
|
---|
See also
-, for, -
Examples
Commpute two different sets of indices to access 4x4 arrays, one for the lower triangular part starting at the main diagonal, and one starting two diagonals further right:
>>> il1 = tril_indices(4)
>>> il2 = tril_indices(4,2)
Here is how they can be used with a sample array: >>> a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]]) >>> a array([[ 1, 2, 3, 4],
[ 5, 6, 7, 8], [ 9, 10, 11, 12], [13, 14, 15, 16]])
Both for indexing: >>> a[il1] array([ 1, 5, 6, 9, 10, 11, 13, 14, 15, 16])
And for assigning values: >>> a[il1] = -1 >>> a array([[-1, 2, 3, 4],
[-1, -1, 7, 8], [-1, -1, -1, 12], [-1, -1, -1, -1]])
These cover almost the whole array (two diagonals right of the main one): >>> a[il2] = -10 >>> a array([[-10, -10, -10, 4],
[-10, -10, -10, -10], [-10, -10, -10, -10], [-10, -10, -10, -10]])
Return the indices for the lower-triangle of an (n,n) array.
See tril_indices() for full details.
Parameters : | n : int
k : int, optional
|
---|
Return the indices for the upper-triangle of an (n,n) array.
Parameters : | n : int
k : int, optional
|
---|
See also
-, for, -
Examples
Commpute two different sets of indices to access 4x4 arrays, one for the upper triangular part starting at the main diagonal, and one starting two diagonals further right:
>>> iu1 = triu_indices(4)
>>> iu2 = triu_indices(4,2)
Here is how they can be used with a sample array: >>> a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]]) >>> a array([[ 1, 2, 3, 4],
[ 5, 6, 7, 8], [ 9, 10, 11, 12], [13, 14, 15, 16]])
Both for indexing: >>> a[iu1] array([ 1, 2, 3, 4, 6, 7, 8, 11, 12, 16])
And for assigning values: >>> a[iu1] = -1 >>> a array([[-1, -1, -1, -1],
[ 5, -1, -1, -1], [ 9, 10, -1, -1], [13, 14, 15, -1]])
These cover almost the whole array (two diagonals right of the main one): >>> a[iu2] = -10 >>> a array([[ -1, -1, -10, -10],
[ 5, -1, -1, -10], [ 9, 10, -1, -1], [ 13, 14, 15, -1]])
Return the indices for the lower-triangle of an (n,n) array.
See triu_indices() for full details.
Parameters : | n : int
k : int, optional
|
---|
Changes consecutive jumps larger than pi to their 2*pi complement.
Pad a time-series with zeros on either side, depending on its length
Parameters : | time_series : n-d array
NFFT : int
|
---|
Returns the z-score of each point of the time series along a given axis of the array time_series.
Parameters : | time_series : ndarray
axis : int, optional
Returns : _______ : zt : ndarray
|
---|