RQA_functions

Features

binscalc()

binscalc(X, n, d, method)

Function to calculate number of bins in a multidimensional histogram using a generalised Freedman–Diaconis rule.

Parameters:

• x (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

Returns:

Bins (array) – number of bins along each dimension

References

• Freedman, David, and Persi Diaconis. “On the histogram as a density estimator: L 2 theory.” Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 57.4 (1981): 453-476.

• Birgé, Lucien, and Yves Rozenholc. “How many bins should be put in a regular histogram.” ESAIM: Probability and Statistics 10 (2006): 24-45.

mutualinfo()

mutualinfo(X, Y, n, d, method='histdd')

Function to calculate mutual information between two time series

Parameters:

• x (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• y (ndarray) – double array of shape (n,d). second time series

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• method (Option between computing the mutual information using:) –

  • multidimensional histogram(“histdd”)

  • average mutual information across dimensions(“avg”)

Returns:

MI (double) – mutual information between time series

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

KNN_MI_vectorized()

KNN_MI_vectorized(X, Y, nearest_neighbor=5, dtype=<class 'numpy.float64'>)

Function to calculate mutual information between two time series using KNN method for datasets that can’t be handled with default binning method. Vectorized version

Parameters:

• X (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• Y (ndarray) – double array of shape (n,d). second time series

• nearest_neighbor (int) – number of nearest neighbour for calculating mutual information, default = 5

Returns:

MI (double) – mutual information between time series

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

• Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 69(6), 066138.

KNN_MI_partial_vectorized()

KNN_MI_partial_vectorized(X, Y, nearest_neighbor=5, dtype=<class 'numpy.float64'>)

Function to calculate mutual information between two time series using KNN method for datasets that can’t be handled with default binning method. Partially version. Vectorized version is faster, however, if size of the time series is large and number of dimensions are much larger, the resulting matrix can’t be stored in the physical memory of the system (RAM) depending on the resource available. In that case this version can be used to use a relatively faster version compared to non-vectorized version

Parameters:

• X (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• Y (ndarray) – double array of shape (n,d). second time series

• nearest_neighbor (int) – number of nearest neighbour for calculating mutual information, default = 5

Returns:

MI (double) – mutual information between time series

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

• Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 69(6), 066138.

KNN_MI_non_vectorized()

KNN_MI_non_vectorized(X, Y, nearest_neighbor=5)

Function to calculate mutual information between two time series using KNN method for datasets that can’t be handled with default binning method. Non-vectorized version. Vectorized version is faster, however, if size of the time series is large and number of dimensions are much larger, the resulting matrix can’t be stored in the physical memory of the system (RAM) depending on the resource available. In that case this version can be used

Parameters:

• X (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• Y (ndarray) – double array of shape (n,d). second time series

• nearest_neighbor (int) – number of nearest neighbour for calculating mutual information, default = 5

Returns:

MI (double) – mutual information between time series

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

• Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 69(6), 066138.

KNN_MI()

KNN_MI(X, Y, nearest_neighbor=5, method='auto', dtype=<class 'numpy.float64'>, memory_limit=4)

Function to calculate mutual information between two time series using KNN method for datasets that can’t be handled with default binning method. Uses vectorised or non-vectorized version depending on whether the required matrix size is less than the specified memory limit

Parameters:

• X (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• Y (ndarray) – double array of shape (n,d). second time series

• nearest_neighbor (int) – number of nearest neighbour for calculating mutual information, default = 5

• method (str) – Specifying options for computing the mutual information using the KNN method. Options are:

  • “auto” : This will check for two additional variables, namely “dtype” (data type) and “memory_limit” (maximum memory allocated for the operation). If, for a given matrix size and data type, the vectorized algorithm fits within the memory limit, it will proceed with the vectorized version. Otherwise, it will compute sequentially.

  • “vectorized” : This will use the vectorized method by default, without checking the memory requirement and the limit specified. This option is faster by default.

  • “sequential” : The algorithm is implemented with for loops instead of vectorization. This could be significantly slower than the vectorized version. However, if resources (RAM/physical memory) are limited and can’t handle huge matrices, this option should be chosen.

  • “partial” : best of both worlds between “vectorized” and “sequential”

• dtype (dtype) – data type, default = np.float64

• memory_limit (double) – memory limit in GiB, default = 4

Returns:

MI (double) – mutual information between time series

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

• Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 69(6), 066138.

timedelayMI()

timedelayMI(u, n, d, tau, method='histdd')

Function to calculate mutual information between a time series and a delayed version of itself

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

Returns:

MI (double) – mutual information between u and u delayed by tau

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

KNN_timedelayMI()

KNN_timedelayMI(u, tau, nearest_neighbor=5, method='auto', dtype=<class 'numpy.float64'>, memory_limit=4)

Function to calculate mutual information between a time series and a delayed version of itself

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• tau (int) – amount of delay

• nearest_neighbor (int) – number of nearest neighbour for calculating mutual information, default = 5

• method (str) – Specifying options for computing the mutual information using the KNN method. Options are:

  • “auto” : This will check for two additional variables, namely “dtype” (data type) and “memory_limit” (maximum memory allocated for the operation). If, for a given matrix size and data type, the vectorized algorithm fits within the memory limit, it will proceed with the vectorized version. Otherwise, it will compute sequentially.

  • “vectorized” : This will use the vectorized method by default, without checking the memory requirement and the limit specified. This option is faster by default.

  • “sequential” : The algorithm is implemented with for loops instead of vectorization. This could be significantly slower than the vectorized version. However, if resources (RAM/physical memory) are limited and can’t handle huge matrices, this option should be chosen.

  • “partial” : best of both worlds between “vectorized” and “sequential”

• dtype (dtype) – data type, default = np.float64

• memory_limit (double) – memory limit in GiB, default = 4

Returns:

MI (double) – mutual information between time series

References

• Shannon, Claude Elwood. “A mathematical theory of communication.” The Bell system technical journal 27.3 (1948): 379-423.

• Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 69(6), 066138.

findtau_default()

findtau_default(u, n, d, grp, mi_method='histdd')

Function to calculate correct delay for estimating embedding dimension based on the first minima of the tau vs mutual information curve

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

Returns:

tau (double) – optimal amount of delay for which mutual information reaches its first minima(and global minima, in case first minima doesn’t exist)

findtau_polynomial()

findtau_polynomial(u, n, d, grp, mi_method='histdd')

Function to calculate correct delay for estimating embedding dimension based on the first minima of the polynomial fit of tau vs mutual information curve

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• grp (str) – group name for saving the TAU vs MI figure

Returns:

tau (double) – optimal amount of delay for which mutual information reaches its first minima(and global minima, in case first minima doesn’t exist)

findtau()

findtau(u, n, d, grp, method='default', mi_method='histdd')

Function to calculate correct delay for estimating embedding dimension based on either the first minima of the tau vs mutual information curve or the polynomial fit of tau vs mutual information curve

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• method (method that should be used to find the first local minima in MI vs tau curve.) – “default” - first minima of the MI vs TAU plot “polynomial”- first minima of the polynomial fit to the MI vs TAU plot

Returns:

tau (double) – optimal amount of delay for which mutual information reaches its first minima(and global minima, in case first minima doesn’t exist)

delayseries()

delayseries(u, n, d, tau, m)

Function to calculate correct delay for estimating embedding dimension based on the first minima of the tau vs mutual information curve

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

• m (int) – number of embedding dimensions

Returns:

s (ndarray) – 3D matrix ov size ((n-(m-1)*tau,m,d)), time delayed embedded signal

References

• Takens, F. (1981). Dynamical systems and turbulence. Warwick, 1980, 366–381.

nearest()

nearest(s, n, d, tau, m)

Function that would give a nearest neighbour map, the output array(nn) stores indices of nearest neighbours for index values of each observations

Parameters:

• s (ndarray) – 3D matrix ov size ((n-(m-1)*tau,m,d)), time delayed embedded signal

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

• m (int) – number of embedding dimensions

Returns:

nn (array) – an array denoting index of nearest neighbour for each observation

References

• Takens, F. (1981). Dynamical systems and turbulence. Warwick, 1980, 366–381.

fnnratio()

fnnratio(u, n, d, m, tau, r, sig)

Function that calculates the ratio of false nearest neighbours

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

• m (int) – number of embedding dimensions

• r (double) – ratio parameter

• sig (double) – standard deviation of the data

Returns:

FNN (double) – ratio of false nearest neighbours, for a given embedding dimension m, when compared an embedding dimension m+1

References

• Kennel, M. B., Brown, R., & Abarbanel, H. D. (1992). Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical review A, 45 (6), 3403.

fnnhitszero()

fnnhitszero(u, n, d, m, tau, sig, delta, Rmin, Rmax, rdiv)

Function that finds the value of r at which the FNN ratio can be effectively considered as zero

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

• m (int) – number of embedding dimensions

• r (double) – ratio parameter

• sig (double) – standard deviation of the data

• delta (double) – the tolerance value(the maximum difference from zero) for a value of FNN ratio to be effectively considered to be zero

• Rmin (double) – minimum value of r from where we would start the parameter search

• Rmax (double) – maximum value of r for defining the upper limit of parameter search

• rdiv (Int) – number of divisions between Rmin and Rmax for parameter search

Returns:

r (double) – value of r at which the value of FNN ratio effectively hits zero

References

• Kantz, H., & Schreiber, T. (2004). Nonlinear time series analysis (Vol. 7). Cambridge university press. section 3.3.1

findm()

findm(u, n, d, tau, sd, delta, Rmin, Rmax, rdiv, bound)

Function that finds the value of m whre the r at which FNN ratio hits zero vs m curve flattens(defined by bound value)

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

• m (int) – number of embedding dimensions

• r (double) – ratio parameter

• sig (double) – standard deviation of the data

• delta (double) – the tolerance value(the maximum difference from zero) for a value of FNN ratio to be effectively considered to be zero

• Rmin (double) – minimum value of r from where we would start the parameter search

• Rmax (double) – maximum value of r for defining the upper limit of parameter search

• rdiv (Int) – number of divisions between Rmin and Rmax for parameter search

• bound (double) – bound value for terminating the parameter serch

Returns:

m (double) – value of embedding dimension

References

• Kantz, H., & Schreiber, T. (2004). Nonlinear time series analysis (Vol. 7). Cambridge university press. section 3.3.1

reccplot()

reccplot(u, n, d, m, tau, eps)

Function that computes the recurrence plot

Parameters:

• u (ndarray) – double array of shape (n,d). Think of it as n points in a d dimensional space

• n (int) – number of samples or observations in the time series

• d (int) – number of measurements or dimensions of the data

• tau (int) – amount of delay

• m (int) – number of embedding dimensions

• eps (double) – radius of the neighbourhood for computing recurrence

Returns:

rplot (ndarray) – recurrence plot

References

• Eckmann, J.-P., Kamphorst, S. O., Ruelle, D., et al. (1995). Recurrence plots of dynamical systems. World Scientific Series on Nonlinear Science Series A, 16, 441–446.

reccrate()

reccrate(rplot, n)

Function that computes the recurrence plot

Parameters:

• rplot (ndarray) – recurrence plot

• n (int) – length of RP

Returns:

reccrate (double) – recurrence rate

References

• Eckmann, J.-P., Kamphorst, S. O., Ruelle, D., et al. (1995). Recurrence plots of dynamical systems. World Scientific Series on Nonlinear Science Series A, 16, 441–446.

findeps()

findeps(u, n, d, m, tau, reqrr, rr_delta, epsmin, epsmax, epsdiv)

Function that computes the recurrence plot

Parameters:

• u (ndarray) – multidimensional time series data

• n (int) – number of observations

• d (int) – number of dimensions

• tau (int) – amount of delay

• m (int) – embedding dimension

• reqrr (doubld) – required recurrence rate specified in the input

• rr_delta (double) – tolerance value for considering a value of recurrence rate to be same as the one that is specified in reqrr

• epsmin (double) – lower bound for the parameter search for epsilon(neighbourhood radius)

• epsmax (double) – upper bound for the parameter search for epsilon(neighbourhood radius)

• epsdiv (double) – number of divisions for the parameter search for epsilon(neighbourhood radius) between epsmin and epsmax

Returns:

eps (double) – epsilon(neighbourhood radius)

References

• Eckmann, J.-P., Kamphorst, S. O., Ruelle, D., et al. (1995). Recurrence plots of dynamical systems. World Scientific Series on Nonlinear Science Series A, 16, 441–446.

vert_hist()

vert_hist(M, n)

Function to compute vertical line distribution(counts of line lengths)

Parameters:

• m (ndarray) – recurrence plot

• n (int) – length of RP

Returns:

nvert (array) – an array containing counts of line lengths

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.

onedhist()

onedhist(M, n)

diaghist()

diaghist(M, n)

Function to compute diagonal line distribution(counts of line lengths)

Parameters:

• m (ndarray) – recurrence plot

• n (int) – length of RP

Returns:

nvert (array) – an array containing counts of line lengths

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.

percentmorethan()

percentmorethan(hst, mini, n)

Function to compute determinism and laminarity from the histogram distribution of lines

The idea is to see what fraction of rcurrent points are part of a linear structure, which can be either vertical or horizontal. The definition of a linear structure on an RP is based on the minimum length(mini) given as an input

Parameters:

• hst (array) – histogram counts of line lengths

• mini (int) – minimum length of consecutive occurances of value 1 in the RP(either vertically or horizontally) that is considered as a line

• n (int) – length of RP

Returns:

nvert (array) – an array containing counts of line lengths

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.

mode()

mode(hst, mini, n)

Function to find mode of the line distributions

Parameters:

• hst (array) – histogram counts of line lengths

• mini (int) – minimum length of consecutive occurances of value 1 in the RP(either vertically or horizontally) that is considered as a line

• n (int) – length of RP

Returns:

p (int) – mode of line length distribution

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.

maxi()

maxi(hst, mini, n)

maximum value in the line length distribution

Parameters:

• hst (array) – histogram counts of line lengths

• mini (int) – minimum length of consecutive occurances of value 1 in the RP(either vertically or horizontally) that is considered as a line

• n (int) – length of RP

Returns:

plmax (int) – max of line length distribution

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.

average()

average(hst, mini, n)

Function to find mean of the line distributions

Parameters:

• hst (array) – histogram counts of line lengths

• mini (int) – minimum length of consecutive occurances of value 1 in the RP(either vertically or horizontally) that is considered as a line

• n (int) – length of RP

Returns:

mu (double) – mean of line length distribution

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.

entropy()

entropy(hst, mini, n)

Function to find entropy of the line distributions

Parameters:

• hst (array) – histogram counts of line lengths

• mini (int) – minimum length of consecutive occurances of value 1 in the RP(either vertically or horizontally) that is considered as a line

• n (int) – length of RP

Returns:

ent (int) – entropy of line length distribution

References

• Webber Jr, C. L., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of applied physiology, 76 (2), 965–973.

• Webber Jr, C. L., & Zbilut, J. P. (2005). Recurrence quantification analysis of nonlinear dynamical systems. Tutorials in contemporary nonlinear methods for the behavioral sciences, 94 (2005), 26–94.

• Marwan, N., Romano, M. C., Thiel, M., & Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Physics reports, 438 (5-6), 237–329.

• Marwan, N., Schinkel, S., & Kurths, J. (2013). Recurrence plots 25 years later—gaining confidence in dynamical transitions. Europhysics Letters, 101 (2), 20007.

• Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., & Kurths, J. (2002). Recurrence- plot-based measures of complexity and their application to heart-rate-variability data. Physical review E, 66 (2), 026702.