From 874bb9fb01dc4137902c35aece2cffe7a300759b Mon Sep 17 00:00:00 2001 From: "Moreno I. Coco" Date: Wed, 3 Dec 2014 00:00:00 +0000 Subject: [PATCH] version 1.0.5 --- DESCRIPTION | 8 +- MD5 | 26 ++-- NAMESPACE | 6 +- R/checkts.R | 54 +++++-- R/crqa.R | 151 ++++++++++++++------ R/optimizeParam.R | 331 +++++++++++++++++++++++-------------------- R/runcrqa.R | 35 +++-- R/wincrqa.R | 90 ++++++++---- man/checkts.Rd | 22 ++- man/crqa.Rd | 52 ++++--- man/optimizeParam.Rd | 36 +++-- man/runcrqa.Rd | 10 +- man/tt.Rd | 2 +- man/wincrqa.Rd | 42 +++--- 14 files changed, 528 insertions(+), 337 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 98a3475..b92db6a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: crqa Type: Package Title: Cross-Recurrence Quantification Analysis for Categorical and Continuous Time-Series -Version: 1.0.4 +Version: 1.0.5 Date: 2014-12-03 Author: Moreno I. Coco and Rick Dale Maintainer: Moreno I. Coco @@ -14,12 +14,12 @@ Description: at the diagonal recurrent points, as well as more in-depth measures of the whole cross-recurrence plot, e.g., recurrence rate. -Depends: R (>= 3.0.0), Matrix, tseriesChaos, fields +Depends: R (>= 3.0.0), Matrix, tseriesChaos, fields, pracma License: GPL (>= 2) Collate: 'CTcrqa.R' 'calcphi.R' 'checkts.R' 'crqa.R' 'drpdfromts.R' 'optimizeParam.R' 'runcrqa.R' 'simts.R' 'spdiags.R' 'takephi.R' 'theiler.R' 'tt.R' 'wincrqa.R' 'windowdrp.R' -Packaged: 2014-03-17 11:30:57 UTC; root +Packaged: 2014-12-05 16:02:34 UTC; root NeedsCompilation: yes Repository: CRAN -Date/Publication: 2014-03-17 12:54:38 +Date/Publication: 2014-12-05 17:20:22 diff --git a/MD5 b/MD5 index 91a8dd8..6ee1bd6 100644 --- a/MD5 +++ b/MD5 @@ -1,36 +1,36 @@ -6d5a3c039a49659684234ddbfca30662 *DESCRIPTION -1cee7a78a67f9d3c8cd68e9ffc4ade52 *NAMESPACE +6ad732bc193dc083039fc23147070c10 *DESCRIPTION +68d8282a548c134580427587d2cf0dd2 *NAMESPACE 28ae7e1d4c946f1c5b80da067e9d2586 *R/CTcrqa.R fd4d3e5903396c9dd00c5d77bb712d2c *R/calcphi.R -b51bdcb8367be865d0aa541d9a4b3649 *R/checkts.R -341931e900fd6808faac508992e58472 *R/crqa.R +0c18c004eda963abd2e0271cc4251b8a *R/checkts.R +15db939b0e514af9d138a38d58c2d495 *R/crqa.R df7954f7b92ca063bb3bd4fff8f61d2d *R/drpdfromts.R -fc7bb8a825747019e0e789157d3eb4b3 *R/optimizeParam.R -dcf5e486656a61ad8de9f92095ff67f2 *R/runcrqa.R +5a3ee97ef5058600237d2a04406b0759 *R/optimizeParam.R +03e91b59adb0dd7aa1c13ea99940c28b *R/runcrqa.R 0264d4c328a7a32ec30714b4a3e16dc8 *R/simts.R ef263b24aca130d3cc80288cd67e1fd4 *R/spdiags.R 47cac4ad60709f2d3a8f9d1b24008fe5 *R/takephi.R 825a40d05a0320a77654ab2475cd6076 *R/theiler.R 5a30c10c022e3bbaa9b9cd235645703d *R/tt.R -3d42f84b035d186565120616e33b9f34 *R/wincrqa.R +02f938031565a9c9fd8295fc161e544b *R/wincrqa.R 54d09d840eb3679de8de57fb308a0b84 *R/windowdrp.R da138fa547a6838629fe7d5b93ba8e90 *data/crqa.RData 9b27c14007b4dfde4894a7a64c1f654d *man/CTcrqa.Rd 5fe54e1726d6996050405a62fdf76285 *man/RDts1.Rd 955a46841cd8af86d1bf7a68537e9d51 *man/RDts2.Rd 96fbb17a2ce0be88599768b0dcea5bc9 *man/calcphi.Rd -5c81c273eaf557d35872ea0919edd910 *man/checkts.Rd +493407e8564b6313bd0d8c4528c10cc5 *man/checkts.Rd 98ae19fd1be11988f8bfe5f49d9c67fb *man/crqa-package.Rd -99f7ac521d1511f667ac9e2d3c3b0456 *man/crqa.Rd +73a17ec5140ba1134be7ce1f30d3ae3d *man/crqa.Rd 212ae1a69c55e5fa7968d72a095be035 *man/drpdfromts.Rd ef2141e99f7f3f2cff924104b4f7fc8e *man/leftmov.Rd -27dbde773468780efc5ad7ea06d1cf44 *man/optimizeParam.Rd +dcc10177c107b376e5eca66be1bad609 *man/optimizeParam.Rd 6548ea97f5c7480527001f75262cba19 *man/rightmov.Rd -068ea7765083d3e35a109cd2cd53cc6d *man/runcrqa.Rd +5afa5d71e31f041413a679f86cc72a71 *man/runcrqa.Rd a02045e8733e9ba46c6d094654cfdb68 *man/simts.Rd 47cbd9f5cd8e35b18bbaf1039e858951 *man/spdiags.Rd a7bbf18e5f6062c94babf3dc8f9c9e1e *man/theiler.Rd -8522ac7eb9dc1e0be3e0066802519f10 *man/tt.Rd -9bac65f30a31c37768464f8e3b01da37 *man/wincrqa.Rd +fafc650cd83a882a5004670ad0598a27 *man/tt.Rd +c03845b1bf61a761d159f45f18d51fed *man/wincrqa.Rd 29dcbbff24f2496573bed92ddffe3810 *man/windowdrp.Rd a9de526d2dd2130a3b2c4eae51b6c365 *src/jspd.f diff --git a/NAMESPACE b/NAMESPACE index c6c4a41..bdcb476 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,7 +15,7 @@ export(CTcrqa, windowdrp) import(tseriesChaos, - fields) + fields, + pracma) -importFrom(Matrix,Matrix) -importFrom(Matrix,t) +importFrom("Matrix", diag, t, Matrix, sparseMatrix, rBind) diff --git a/R/checkts.R b/R/checkts.R index c89e554..fccf4c9 100644 --- a/R/checkts.R +++ b/R/checkts.R @@ -1,15 +1,17 @@ ## arguments: ts1, ts2, the two time-series, ## datatype: (numerical, categorical) -## +## pad: whether the two series of to be chopped or padded -checkts <- function(ts1, ts2, datatype, thrshd){ +.packageName <- 'crqa' + +checkts <- function(ts1, ts2, datatype, thrshd, pad){ ## take as input the two sequences and normalize them by a pre-defined ## threshold constant. ln1 = length(ts1) ln2 = length(ts2) - + mxl = max(c(ln1, ln2)) if (datatype == "numeric"){ ## just make sure that the vec are numeric ts1 = as.numeric(as.matrix( ts1) ) @@ -28,17 +30,41 @@ checkts <- function(ts1, ts2, datatype, thrshd){ } if ( ln1 != ln2 ){ ## timeseries can differ in length - df = abs(ln1 - ln2) + dfs = abs(ln1 - ln2) ## TODO if length of time-series is different there is a threshold ## find optimal solution to accepted threshold - if (df <= thrshd){ ## threshold to adjust | remove pair - - if (ln2 > ln1){ ts2 = ts2[1:(length(ts2)-df)] } ## remove final time-points to series of equal length - if (ln1 > ln2){ ts1 = ts1[1:(length(ts1)-df)] } - - } + if (dfs <= thrshd){ ## threshold to adjust | remove pair + + ## remove final time-points to series of equal length + + if (pad == TRUE){ + if (datatype == "numeric"){ + ## we pad with the mean value + + if (ln2 > ln1){ ts1 = c(ts1, rep(mean(c(ts1,ts2), na.rm = TRUE), dfs)) } + if (ln1 > ln2){ ts2 = c(ts2, rep(mean(c(ts1,ts2), na.rm = TRUE), dfs)) } + + } + + if (datatype == "categorical"){ + ## we pad with a constant that it is not used + ## for recoding + + if (ln2 > ln1){ ts1 = c(ts1, rep(mxl + 1, dfs)) } + if (ln1 > ln2){ ts2 = c(ts2, rep(mxl + 1, dfs)) } + + } + + + } else { + ## we just chop off the series + if (ln2 > ln1){ ts2 = ts2[1:(length(ts2)-dfs)] } + if (ln1 > ln2){ ts1 = ts1[1:(length(ts1)-dfs)] } + } + + } } if ( length(ts1) == length(ts2) ){ @@ -47,7 +73,7 @@ checkts <- function(ts1, ts2, datatype, thrshd){ } else { - return (list(df, FALSE) )## how many units are the sequences differing. + return (list(dfs, FALSE) )## how many units are the sequences differing. } @@ -65,8 +91,10 @@ numerify <- function( ts1, ts2 ){ ids = seq(1, length(objects), 1) for ( o in 1:length(objects) ){ - nwts1[which(ts1 == objects[o])] = ids[o] - nwts2 [which(ts2 == objects[o])] = ids[o] + indts1 = which(ts1 == objects[o]) + if (length(indts1) > 0) nwts1[indts1] = ids[o] + indts2 = which(ts2 == objects[o]) + if (length(indts2) > 0) nwts2[indts2] = ids[o] } return( list( nwts1, nwts2 ) ) diff --git a/R/crqa.R b/R/crqa.R index fb9ad8e..1daafd3 100644 --- a/R/crqa.R +++ b/R/crqa.R @@ -1,6 +1,6 @@ ## written in R by Moreno I. Coco, 2013, (moreno.cocoi@gmail.com) ## crqa, adapted from a Matlab code developed at -## summer school Nonlinear Methods for Psychological Science +## summer school of: Nonlinear Methods for Psychological Science ## organized by the University of Cincinnati, 2012 ## next time i check change sd(matrix) with sapply(matrix, sd) @@ -10,72 +10,104 @@ ## embed = the embedding dimension, i.e., the lag intervals ## rescale = the normalization for the distance; ## if 1 (Mean Distance); if 2 (Max Distance) -## radius = the maximum distance to be calculated (set it very -## small, if the series are categorical in nature +## radius = distance to accept two points as recurrent (set it very +## small, if the series are categorical in nature) ## normalize = rescale factor for source data; ## if 1 (Unit interval); if 2 (z-score) ## mindiagline = set a minimum diagonal line length ## mindiagline = set a minimum vertical line length -# delay = 1; embed = 1; rescale = 1; radius = 0.001; -# normalize = 0; mindiagline = 2; minvertline = 2; -# whiteline = FALSE # - flag to compute or not white vertical lines -# in the recurrence plot. Note, white lines are not +## whiteline = FALSE # - flag to compute or not white vertical lines +## in the recurrence plot. Note, white lines are not ## yet used to derive any particular measure -# recpt = FALSE # - flag to indicate whether the input ts1 is already +## recpt = FALSE # - flag to indicate whether the input ts1 is already ## a recurrence plot + ## tw = the size of the Theiler Window, the default is 0 +## side = a string indicating whether the recurrence measures +## should be calculated in the "upper" triangle of the matrix +## "lower" triangle of the matrix, on the "whole" matrix + +## checkl = a list with four arguments: +## $do = TRUE|FALSE; normalize (or not) the length of ts +## if $do = TRUE, then the arguments of checkts() needs to be passed. +## $datatype = (numerical, categorical) - nature of ts +## $thrshd: number of timepoints we accept two ts to be different +## $pad: whether the two series have to be padded (TRUE) or chopped + +## try below + ## ts1 = c(0,0,1,1,0,0) ## ts2 = c(2,2,1,1,2,2) -## ts1 = c(0, 0, 1, 1, 0, 0, 2, 2, 1, 1) -## ts2 = c(2, 2, 1, 1, 2, 2, 1, 1, 0, 0) +# ts1 = c(0, 0, 1, 1, 0, 0, 2, 2, 1, 1) +# ts2 = c(1,1, 2, 2, 0, 0, 1, 2) +#delay = 1; embed = 1; rescale = 1; radius = 0.001; +#normalize = 0; mindiagline = 2; minvertline = 2; +#tw = 0; whiteline = FALSE; recpt = FALSE; side = "both" +#checkl = list(do = TRUE, thrshd = 3, datatype = "categorical", +# pad = TRUE) -packageName <- 'crqa' +#crqa(ts2, ts1, delay, embed, rescale, radius, normalize, mindiagline, minvertline, tw, whiteline, recpt, side, checkl) + +.packageName <- 'crqa' crqa <- function(ts1, ts2, delay, embed, rescale, radius, normalize, mindiagline, minvertline, - tw, whiteline, recpt){ + tw = 0, whiteline = F, recpt = F, side = "both", + checkl = list(do = F)){ - if( missing(tw) ){ tw = 0} ## default for Theiler window + ## passing a few default above + # if( missing(tw) ){ tw = 0} -> not working v11 = v21 = NULL ## stupid initializations to please CRAN - ## require("fields") ## to compute the Euclidean distance matrix - ## require("Matrix") ## to manipulate sparse matrices + # require("fields") ## to compute the Euclidean distance matrix + # require("Matrix") ## to manipulate sparse matrices ## check if the input is a recurrence plot if (recpt == FALSE){ - ts1 = as.vector(as.matrix(ts1)) ## make sure data is a vector ts2 = as.vector(as.matrix(ts2)) if (is.matrix(ts1)){ stop("Your data must consist of a single column of data.")} - if (is.matrix(ts2)){ stop("Your data must consist of a single column of data.")} - - - ##chop of sequences if they are of different lengths - - if (length(ts1) != length(ts2)){ - shortest = min(c(length(ts1), length(ts2)) ); - ts1 = ts1[1:shortest]; - ts2 = ts2[1:shortest]; + if (is.matrix(ts2)){ stop("Your data must consist of a single column of data.")} + + ## check the length of the sequences and decide if they have + ## to be normalized to the same length. + + if (checkl$do == TRUE){ + + tsnorm = checkts(ts2, ts1, checkl$datatype, + checkl$thrshd, checkl$pad) + + if (tsnorm[[2]] == FALSE){ + stop("Time-series difference longer than threshold. Increase threshold, or set checkl$do = FALSE avoiding normalization of ts") + } else { + ts1 = tsnorm[[1]][,1] + ts2 = tsnorm[[1]][,2] + } + } + + - ##rescale the data if really necessary + ##rescale the data if really necessary if (normalize > 0){ switch (normalize, - {1 + {1 + ## unit-interval ts1 = (ts1 - min(ts1)); ts1 = ts1 / max(ts1); ts2 = (ts2 - min(ts2)); ts2 = ts2 / max(ts2);}, - {2 - ts1 = (ts1 - mean(ts1))/sd(ts1) # zscore(ts1, na.rm = TRUE, robust = FALSE); ## using R.basics - ts2 = (ts2 - mean(ts2))/sd(ts2) # zscore(ts2, na.rm = TRUE, robust = FALSE) + {2 + ## z-score + ts1 = (ts1 - mean(ts1))/sd(ts1) + ts2 = (ts2 - mean(ts2))/sd(ts2) } ) } @@ -127,7 +159,11 @@ crqa <- function(ts1, ts2, delay, embed, rescale, ## Compute euclidean distance matrix - vlength = length(v11) ## just to have the length of matrix saved + v1l = length(v11) + v2l = length(v21) + + + ## just to have the length of matrix saved dm = rdist(dimts1,dimts2); ## Find indeces of the distance matrix that fall @@ -149,13 +185,18 @@ crqa <- function(ts1, ts2, delay, embed, rescale, } else { dmrescale = dm } ## Compute recurrence matrix - ind = which(dmrescale <= radius, arr.ind=T); + ind = which(dmrescale <= radius, arr.ind = TRUE); r = ind[,1]; c = ind[,2] } else { ## take as input an RP directly - vlength = nrow(ts1) - ind = which(ts1 > 0, arr.ind = T) + ## as usual R needs fiddly code to make sure about identify of data + ts1 = matrix(as.logical(ts1), ncol = ncol(ts1)) + v1l = nrow(ts1); v2l = ncol(ts1) + + ## matrix needs to be logical + ind = which(ts1 > 0, arr.ind = TRUE) + ## just a trick to reduce the number of lines ## of the code r = ind[,1]; c = ind[,2] @@ -163,7 +204,8 @@ crqa <- function(ts1, ts2, delay, embed, rescale, } if (length(r) != 0 & length(c) != 0){ ##avoid cases with no recurrence - S = sparseMatrix(r, c, dims = c(vlength, vlength)) + + S = sparseMatrix(r, c, dims = c(v1l, v2l)) ## this is the recurrent plot ## transpose it to make identical to Marwan S = t(S) @@ -172,13 +214,32 @@ crqa <- function(ts1, ts2, delay, embed, rescale, ## Marwan blanks out the recurrence along the diag S = theiler(S, tw) + if (side == "upper"){ + ## if only the upper side is of interest + ## it blanks out the lowest part + S = as.matrix(S) + S[lower.tri(S, diag = TRUE)] = 0 + S = Matrix(S, sparse = TRUE) + } + + if (side == "lower"){ + ## viceversa + S = as.matrix(S) + S[upper.tri(S, diag = TRUE)] = 0 + S = Matrix(S, sparse = TRUE) + } + + if (side == "both"){ + ## just keep it as is. + S = S} + spdiagonalize = spdiags(S) ## spdiags should have decent speed B = spdiagonalize$B ##calculate percentage recurrence by taking all non-zeros numrecurs = length(which(B == TRUE)); - percentrecurs = (numrecurs/((vlength^2)))*100; + percentrecurs = (numrecurs/((v1l*v2l)))*100; #################################################################### #################################################################### @@ -263,16 +324,16 @@ crqa <- function(ts1, ts2, delay, embed, rescale, lam = 0; TT = 0; RP = NA } - results = list(rec = percentrecurs, det = pdeter, - nrline = numdiaglines, maxline = maxline, - meanline = meanline, entropy = entropy, - relEntropy = relEntropy, - lam = lam, TT = TT, RP = S) + results = list(RR = percentrecurs, DET = pdeter, + NRLINE = numdiaglines, maxL = maxline, + L = meanline, ENTR = entropy, + rENTR = relEntropy, + LAM = lam, TT = TT, RP = S) - } else { print (paste ("No recurrence found") ) - results = list(rec = 0, det = NA, nrline = 0, maxline = 0, - meanline = 0, entropy = NA, relEntropy = NA, - lam = NA, TT = NA, RP = NA)} + } else { # print (paste ("No recurrence found") ) + results = list(RR = 0, DET = NA, NRLINE = 0, + maxL = 0, L = 0, ENTR = NA, rENTR = NA, + LAM = NA, TT = NA, RP = NA)} return (results) diff --git a/R/optimizeParam.R b/R/optimizeParam.R index 8283168..361c3df 100644 --- a/R/optimizeParam.R +++ b/R/optimizeParam.R @@ -10,134 +10,97 @@ # TODO: probably, for windows the selection of size/embedding # has to be done at the same time. - -## source("./ami.R") ## average mutual information ## between 2 series ## source("./crqa.R") +# lgM = 20 ## maximum lag to check for average mutual information -## lgM = 20 ## maximum lag to check for average mutual information -## steps = seq(1, 10, 1) ##how many points we should look -## ahead to define a local minimum - -## cut.del = seq(1, 40,1) ## cut off of delay is set -## radiusspan = 2 ## increasing radius step means smaller steps -## radiussample = 10 ## number of radius points within the steps +#radiusspan = 2 ## increasing radius step means smaller steps +#radiussample = 10 ## number of radius points within the steps ## to be sampled +#par = list(lgM = 20, radiusspan = 100, radiussample = 40, +# normalize = 0, rescale = 1, mindiagline = 2, steps = seq(1,6,1), +# minvertline = 2, tw = 0, whiteline = FALSE, fnnpercent = 10, +# recpt = FALSE) + +## min.rec = minimal recurrence value +## max.rec = maximal recurrence value + +#ts1 = runif(100) +#ts2 = runif(100) + .packageName <- 'crqa' -optimizeParam <- function(ts1, ts2, par){ +optimizeParam <- function(ts1, ts2, par, min.rec=2, max.rec=5){ + rescale = normalize = mindiagline = minvertline = - lgM = steps = cut.del = recpt = whiteline = tw = - radiusspan = radiussample = NULL - ## initialize attached variables - - for (v in 1:length(par)) assign(names(par)[v], par[[v]]) ## assign parameters - - if (radiusspan <= 1){ - stop("Radius span too small, please choose a larger value") - } - # require(tseriesChaos) - - ## check that the input is vector and numeric + lgM = fnnpercent = recpt = whiteline = tw = + steps = radiusspan = radiussample = NULL + + ## require(tseriesChaos) + + ## Initialize attached variables + ## Check that the input is vector and numeric if (is.vector(ts1) != TRUE){ ts1 = as.numeric(as.matrix(ts1) ) } if (is.vector(ts2) != TRUE){ ts2 = as.numeric(as.matrix(ts2) ) } - mi1 = as.numeric( mutual(ts1, lag.max = lgM, plot = FALSE)) - mi2 = as.numeric( mutual(ts2, lag.max = lgM, plot = FALSE)) + for (v in 1:length(par)) assign(names(par)[v], par[[v]]) + ## Assign parameters - ## iterate over a range of possible steps to find out the optimal one + if (radiusspan <= 1){ + stop("Radius span too small, please choose a larger value") + } + + ## A: Choose a delay that will accommodate both ts. + ## For example if one has a considerably longer delay + ## indicated than another, you may need to select the longer + ## delay of the two because that ensures that new information + ## is gained for both by using the higher delay. + ## If delays are close to each other, you may want to use a + ## delay somewhere in between the two + + mi = as.numeric( mutual(rbind(ts1, ts2, deparse.level = 0), + lag.max = lgM, plot = FALSE)) + + lag = vector() + s = 1 + + if (length(mi) < length(steps)){ + stop("Decrease number of steps") + } - lag1 = lag2 = vector() for (s in steps){ # print(s) for (l in 1:(lgM - s)){ - if (mi1[l] < mi1[l + s]){ - lg1 = which(mi1 == mi1[l]) - lag1 = c(lag1, lg1) + if (mi[l] < mi[l + s]){ + lg = which(mi == mi[l]) + lag = c(lag, lg) break } } - for (l in 1:(lgM - s)){ - if (mi2[l] < mi2[l + s]){ - lg2 = which(mi2 == mi2[l]) - lag2 = c(lag2, lg2) - break } - } } - if (length(lag1) == 0 | length(lag2) == 0){ + if (length(lag) == 0){ stop("Please try varying maximum lag (and/or step size): minimal mutual information was not reached") } - ## take unique lags, i.e., avoid repeat - lag1 = unique(lag1); lag2 = unique(lag2) - - - ## for each time-series take the lag where MI is minimal - - ix.lg1 = which( mi1[lag1] == min(mi1[lag1]) ) - if (length(ix.lg1) > 1){ - ix.lg1 = sample(ix.lg1, 1) ## if there are more than one minimum - ## sample one - lg1 = lag1[ix.lg1]} - - ix.lg2 = which( mi2[lag2] == min(mi2[lag2]) ) - if (length(ix.lg2) > 1){ - ix.lg2 = sample(ix.lg2, 1) - lg2 = lag2[ix.lg2]} - - -# mi = ami(ts1, ts2, lag = 1)[[1]] ## not sure about this function -# ## perhaps look at it later - -# Choose a delay that will accommodate both ts. -# For example if one has a considerably longer delay -# indicated than another, you may need to select the longer -# delay of the two because that ensures that new information -# is gained for both by using the higher delay. -# If delays are close to each other, you may want to use a -# delay somewhere in between the two. - - delay.in = vector() - - for (delay in cut.del){ - ## print(delay) - - if (lg1 > lg2 & abs(lg1 - lg2) > delay){ - ## print( c(delay, lg1) ) - delay.in = c(delay.in, lg1) - } - - else if (lg2 > lg1 & abs(lg2 - lg1) > delay){ - ## print( c(delay, lg2) ) - delay.in = c(delay.in, lg2) - } - - else { - delay.in = c(delay.in, - round( mean(c(lg1, lg2))) ) - } - - ## print(delay.in) - - } - - ## take the delay that more frequently is found - ## to be optimal - - del.in = as.data.frame( table(delay.in) ) + del.in = as.data.frame( table(lag) ) del = as.numeric( as.vector( del.in[which(del.in[,2] == max(del.in[,2])),1])) - if (length(del) > 1){ del = del[length(del)] } ## take longer delay - - - ## Determine embedding dimension: %FNN function bottoms out + if (length(del) > 1){ del = del[length(del)] } + + del = del-1 + ## to account for the 0 delay + + ## another solution would be to use the global minimum + ## but it could result into longer delays being selected + + ## B: Determine embedding dimension: %FNN function bottoms out ## (i.e., where it buys you nothing to add more dimensions). ## If the embedding dimension indicated for the two files is ## different, you’ll want to go with the higher embedding @@ -148,28 +111,66 @@ optimizeParam <- function(ts1, ts2, par){ ## high embedding dimensions. ## TODO: bring out this parameters to optimize them also -# m maximum embedding dimension = 20 -# d delay parameter = get it from previous optimization -# t Theiler window = 0 ## in cross-recurrence LOC is important -# rt escape factor = leave it default -# eps neighborhood diameter = leave default - + ## m maximum embedding dimension = 20 + ## d delay parameter = get it from previous optimization + ## t Theiler window = 0 + ## in cross-recurrence LOC is important + ## rt escape factor = leave it default + ## eps neighborhood diameter = leave default + embdts1 = false.nearest(ts1, m = 20, d = del, t = 0, rt = 10, eps = sd(ts1)/10) - embmints1 = as.numeric( which(embdts1[1,] == min(embdts1[1,], - na.rm = TRUE) ) ) - + + ## get a percentage of reduction of false neighbours + ## based on the first dimension + + fnnfraction1 = embdts1[1,] + fnnfraction1 = fnnfraction1[which(is.na(fnnfraction1) == FALSE)] + emdthd1 = fnnfraction1[1]/fnnpercent + + emdix1 = which(diff(fnnfraction1) < -emdthd1) + + if (length(emdix1) == 1){ + emdmints1 = as.numeric(emdix1) + 1 + } else if (length(emdix1) > 1){ + ## there is no gain when embedding + emdmints1 = as.numeric(tail(emdix1,1) + 1) + } else { + emdmints1 = 1 + } + + ## to adjust for the diff + + ## alternative method, just get the minimum + ## at the moment the method is commented. + # embmints1 = as.numeric( which(fnnfraction1 == min(fnnfraction1))) + embdts2 = false.nearest(ts2, m = 20, d = del, t = 0, rt=10, eps=sd(ts2)/10) - embmints2 = as.numeric( which(embdts2[1,] == min(embdts2[1,], na.rm = TRUE) ) ) + + fnnfraction2 = embdts2[1,] + fnnfraction2 = fnnfraction2[which(is.na(fnnfraction2) == FALSE)] + emdthd2 = fnnfraction2[1]/fnnpercent + + emdix2 = which(diff(fnnfraction2) < -emdthd2) + + if (length(emdix2) == 1){ + emdmints2 = as.numeric(emdix2) + 1 + } else if (length(emdix2) > 1){ + emdmints2 = as.numeric(tail(emdix2,1) + 1) + } else { ## there is no gain when embedding + emdmints2 = 1 + } + + # embmints2 = as.numeric( which(fnnfraction2 == min(fnnfraction2))) - if ( length(embmints1) > 1){ embmints1 = embmints1[1] } - if ( length(embmints2) > 1){ embmints2 = embmints2[1] } + if ( length(emdmints1) > 1){ emdmints1 = emdmints1[1] } + if ( length(emdmints2) > 1){ emdmints2 = emdmints2[1] } - embdim = max( c(embmints1, embmints2) ) + embdim = max( c(emdmints1, emdmints2) ) - ## take the optimal parameters from above and then check - ## set the radius + ## C: Use the optimal parameters of embedding and delay + ## from above and then explore an optimal radius ## first we need to see in the rescaled matrix what is ## the range of values, so that the radius interval @@ -191,57 +192,85 @@ optimizeParam <- function(ts1, ts2, par){ dmrescale = (dm/rescaledist)*100} ) } else { dmrescale = dm } - - ## take the sd error as a unit for the measure - combo = c(ts1,ts2) + + + combo = c(ts1, ts2) sdun = sd(dmrescale) - mnun = median(dmrescale) ## the distance that gives us RR 50% - - ## sequence of radius from max to 0 in steps - ## relative to SD + mnun = median(dmrescale)*2 + ## Multiplier to increase range of candidate radii + radi = seq(mnun, 0, -(sdun/radiusspan)) - - ## we take only the lower half, where it is certain that - ## radius will produce RR < 25% radi = radi[(length(radi)/2):length(radi)] - - kpt = ceiling(length(radi)/radiussample) rsamples = sample(1:kpt, 1) - - syssamp = seq(rsamples, rsamples + kpt*(radiussample-1), kpt) + + syssamp = seq(rsamples, rsamples + kpt * (radiussample - 1), kpt) syssamp = syssamp[syssamp <= length(radi)] - radi = radi[syssamp] - - ## delay and embed dimension to pass - delay = del; embed = embdim; - delay = 1; embed = 1 - + delay = del + embed = embdim optrad = vector() - # count = 0 - for (r in radi){ - # count = count + 1 - ## print( paste( "Checking radius:", r) ) - radius = r - res = crqa(ts1, ts2, delay, embed, - rescale, radius, normalize, mindiagline, - minvertline, tw, whiteline, recpt) + + end.flag <- 0 + + while(end.flag==0){ + ## location with largest radius + hi.loc<-1 + ## location with smallest radius + lo.loc<-length(radi) + + #print(c("hi.loc",hi.loc)) + #print(c("lo.loc",lo.loc)) + ## take middle location as the current readius to be tested - # print(res$rec) + curr.loc<-round(length(radi)/2) + radi[curr.loc]->r + radi + + #print(c("r",r)) + #print(c("curr.loc",curr.loc)) + #print(c("embed",embed)) + #print(c("delay",delay)) - if (res$rec >= 2 & res$rec <= 5){ - optrad = r - break} + res = crqa(ts1, ts2, delay, embed, rescale, r, normalize, + mindiagline, minvertline, tw, whiteline, recpt) + #print(c("recurr",res$RR)) + #print("#######") + #print("#######") - } - - if (length(optrad) == 0){ - optrad = NA} + if (res$RR >=min.rec & res$RR <= max.rec) { + optrad = r + #print(c("radius",optrad)) + #print(c("embed",embed)) + #print(c("delay",delay)) + end.flag<-1 + } + else { + + if (res$RRlo.loc + } + if (res$RR>max.rec){ + ## if rec greater than max rec, curr.loc becomes new hi + curr.loc->hi.loc + } + if((lo.loc-hi.loc)<2){ + ## if less than 2 radi locs remaining, no optimal radius + end.flag<-1 + warning("Optimal Radius Not found: try again choosing a wider radius span and larger sample size") + } + } + ## replace radi vector with remaining unsearched vector + radi<-radi[hi.loc:lo.loc] + } ## end while - if(is.na(optrad)){ - warning("Optimal radius not found: try again choosing a wider radius span and larger sample size")} + if (length(optrad) == 0) { + optrad = NA + } - return ( list(radius = optrad, emddim = embdim, delay = del) ) + if (!is.na(optrad)){ + return(list(radius = optrad, emddim = embdim, delay = del)) + } } diff --git a/R/runcrqa.R b/R/runcrqa.R index 3ffc5ec..34d4010 100644 --- a/R/runcrqa.R +++ b/R/runcrqa.R @@ -21,7 +21,8 @@ ## cons: it does not return the full spectrum of measures ## from the recurrence plots -## type = 2; a recurrence matrix is calculated and several measures of recurrence: meanline, maxline, det, entropy ... are extracted from the matrix +## type = 2; a recurrence matrix is calculated and several measures of +## recurrence: meanline, maxline, det, entropy ... are extracted from the matrix ## again recurrence can be computed both on the whole profile ## and based on the window. @@ -58,25 +59,30 @@ runcrqa <- function(ts1, ts2, par){ datatype = thrshd = type = method = ws = radius = windowsize = - lagwidth = delay = rescale = normalize = mindiagline = minvertline = - tw = whiteline = recpt = NULL + windowstep = delay = rescale = normalize = mindiagline = + minvertline = lagwidth = tw = whiteline = recpt = pad = NULL ## stupid initialization to please CRAN - for (v in 1:length(par)) assign(names(par)[v], par[[v]]) ## assign parameters + for (v in 1:length(par)) assign(names(par)[v], par[[v]]) + ## assign parameters - tryCatch({ ## set up a tryCatch to detach parameters values if errors occur + tryCatch({ + ## set up a tryCatch to detach parameters values if errors occur - res = checkts(ts1, ts2, datatype, thrshd) ## first check that sequences - ## have the same length + res = checkts(ts1, ts2, datatype, thrshd, pad) + ## first check that sequences + ## have the same length if ( res[[2]] == TRUE ){ tsnorm = res[[1]] ts1 = tsnorm[,1]; ts2 = tsnorm[,2] - switch(type, # Ways of calculating recurrence + switch(type, + ## Ways of calculating recurrence - {1 ## Quick Recurrence Profile (only diagonal) + {1 + ## Quick Recurrence Profile (only diagonal) if (method == "profile"){ @@ -100,16 +106,17 @@ runcrqa <- function(ts1, ts2, par){ if (method == "profile"){ res = crqa(ts1, ts2, delay, embed, - rescale, radius, normalize, mindiagline, minvertline, - tw, whiteline, recpt) + rescale, radius, normalize, mindiagline, + minvertline, tw, whiteline, recpt) } if (method == "window"){ - res = wincrqa(ts1, ts2, step, windowsize, lagwidth, delay, - embed, rescale, radius, normalize, mindiagline, minvertline, - tw, whiteline, recpt) + res = wincrqa(ts1, ts2, windowstep, + windowsize, delay, embed, rescale, + radius, normalize, mindiagline, minvertline, + tw, whiteline, trend = F) } diff --git a/R/wincrqa.R b/R/wincrqa.R index 8d9f5cb..934aaba 100644 --- a/R/wincrqa.R +++ b/R/wincrqa.R @@ -1,55 +1,85 @@ -## GNU License, written by Moreno I. Coco + +## written by Moreno I. Coco 2013 (moreno.cocoi@gmail.com) ## original Matlab code by Rick Dale ## calculate recurrence over different sized windows ## arguments: step = interval considered on the serie; -## window_size = the size of the window wherin crqa is runned. -## lag_profile_width = lags within the window -## opt = list( type = 1 - do, or 0 - don't do crqa on windows; -## if type 1; opt should also specify parameters to crqa -## (delay, embed, rescale, radius, normalize, minline) -## -# step = 10; -# windowsize = 100; -# lagwidth = 20; -# type = 1; delay = 1; embed = 1; rescale = 1; radius = 0.00001; normalize = 0; minline = 2 +## windowstep = the step of the window. +## windowsize = the size of the window +## +## +## other arguments to pass are the same as crqa + +# windowsize = 200; windowstep = 50; +# type = 1; delay = 1; embed = 1; +# rescale = 1; radius = 0.0001; +# normalize = 0; minline = 2 # source("crqa.R") +# tS = simts(0.25, 0.05, 0.2, 0.2, 0.25, 1000) +# ts1 = tS[1,]; ts2 = tS[2,] + .packageName <- 'crqa' -wincrqa <- function(x, y, step, windowsize, lagwidth, delay, embed, -rescale, radius, normalize, mindiagline, minvertline, tw, whiteline, recpt){ +wincrqa <- function(ts1, ts2, windowstep, windowsize, delay, embed, + rescale, radius, normalize, mindiagline, + minvertline, tw = 0, whiteline = F, trend = FALSE){ - recpt = FALSE - ## we do not expect as input a windowed recurrent plot - x = as.vector(as.matrix(x)); y = as.vector(as.matrix(y)) - points = seq(1, (length(x) - (windowsize)-1), step) - - crawin = vector() - tsp = 0 ## set a counter with the win at which rec was found - + ## we do not expect as input a recurrent plot + ts1 = as.vector(as.matrix(ts1)); ts2 = as.vector(as.matrix(ts2)) + points = seq(1, (length(ts1) - (windowsize)-1), windowstep) + TREND = crawin = vector() + tsp = 0 ## set a counter with the win at which rec was found + + i = 1 for (i in points){ tsp = tsp +1 - xwin = x[i:(i+windowsize)]; - ywin = y[i:(i+windowsize)]; + ts1win = ts1[i:(i+windowsize)]; + ts2win = ts2[i:(i+windowsize)]; - ans = crqa(xwin, ywin, delay, embed, rescale, + ans = crqa(ts1win, ts2win, delay, embed, rescale, radius, normalize, mindiagline, minvertline, tw, - whiteline, recpt) + whiteline) - + RP = ans$RP + if (length(RP) == 1) RP = vector() ## a trick for cases + ## with empty recurrence plot ans = as.numeric( unlist(ans[1:9]) ) - ## for window recurrent do not save the recurrence plot + ## if trend needs to be calculated do it here + + if (trend == TRUE){ + + if (length(RP) > 0){ + + NX = ncol(RP) + T = vector("numeric", length = NX-1) + + k = 1 + for (k in 1:(NX-1)){ + T[k] = length( + which(diag(RP, k) != F)) / (NX-k)*100; + } + + Ntau = NX - 1 - round(0.1*NX); + + ## last 10% of the RP will be skipped + p = polyfit(2:(Ntau+1),T[1:Ntau], 1) # slope + TREND = c(TREND, 1000 * p[1]) + ## Webber's definition includes factor 1000 + + } else { TREND = NA} + } else { + TREND = c(TREND, NA) + } crawin = rbind(crawin, c(ans, tsp), deparse.level = 0) } - - - return(crawin) + + return(list(crqwin = crawin, TREND = TREND)) } diff --git a/man/checkts.Rd b/man/checkts.Rd index 24da26f..08e120b 100644 --- a/man/checkts.Rd +++ b/man/checkts.Rd @@ -6,11 +6,10 @@ whether they have to be discarded given a specified threshold. If the difference between series is smaller than the threshold, the longest series between - the two is trimmed from the tail to be the same length - as the shortest series.} + the two is either padded, or trimmed, to be the same length.} \usage{ -checkts(ts1, ts2, datatype, thrshd) +checkts(ts1, ts2, datatype, thrshd, pad) } \arguments{ @@ -22,6 +21,13 @@ checkts(ts1, ts2, datatype, thrshd) two-series used to accept them as valid. Series with a difference in length bigger than threshold are rejected as invalid.} + \item{pad}{A logical flag indicating whether, + in case the two-series differ by length, + have to be padded (TRUE), with the shorter being + extended. Or not, in which case, the longest + series is trimmed to be of the same length + as the shortest.} + } \details{ This function strictly applies when the two time-series @@ -32,13 +38,16 @@ checkts(ts1, ts2, datatype, thrshd) to discriminate small versus large differences. The value to assign to the threshold should be estimated by looking at the difference distribution observed - in the dataset. + in the dataset. + Note, when the two series are of continous numerical + values, the padding value is obtained as the mean + between the two series. } \value{If the difference between the two-series is smaller than the threshold, it returns a list with two arguments: the first [[1]] is a matrix with two columns, i.e., the - two time-series trimmed, and number of rows equal to the + two time-series padded or trimmed, and number of rows equal to the shortest sequence. The second [[2]] is a boolean flag with value TRUE, indicating that the difference was shorter than the threshold. If instead the difference @@ -60,8 +69,9 @@ checkts(ts1, ts2, datatype, thrshd) ts1 = seq(1,30,1); ts2 = seq(1,25,1) datatype = "continuous" threshold = 6 ## threshold is larger than difference +pad = FALSE -res = checkts(ts1, ts2, datatype, threshold) +res = checkts(ts1, ts2, datatype, threshold, pad) print(res) threshold = 4 ## threshold is smaller than difference diff --git a/man/crqa.Rd b/man/crqa.Rd index 483da87..077deca 100644 --- a/man/crqa.Rd +++ b/man/crqa.Rd @@ -14,7 +14,7 @@ extracted (explained below).} \usage{ crqa(ts1, ts2, delay, embed, rescale, radius, normalize, -mindiagline, minvertline, tw, whiteline, recpt) +mindiagline, minvertline, tw, whiteline, recpt, side, checkl) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -43,8 +43,17 @@ mindiagline, minvertline, tw, whiteline, recpt) or not (FALSE) empty vertical lines.} \item{recpt}{A logical flag indicating whether measures of cross-recurrence are calculated directly - from a recurrent plot (TRUE) or not (FALSE). - } + from a recurrent plot (TRUE) or not (FALSE).} + \item{side}{A string indicating whether recurrence measures + should be calculated in the 'upper' triangle of the RP + 'lower' triangle of the matrix, or 'both'. + LOC is automatically excluded for 'upper' and 'lower'.} + \item{checkl}{A list with four arguments: + do = TRUE|FALSE; normalize (or not) the length of ts + if do == TRUE, then the arguments of checkts() needs to be passed. + datatype = (numerical, categorical) - nature of ts + thrshd: number of timepoints we tollerate difference between ts + pad: the two series have to be padded (TRUE) or chopped (FALSE)} } \details{ @@ -59,22 +68,22 @@ mindiagline, minvertline, tw, whiteline, recpt) plot. Otherwise, the values for the output arguments will be either 0 or NA. - \item{rec}{The percentage of recurrent points falling within + \item{RR}{The percentage of recurrent points falling within the specified radius (range between 0 and 100)} - \item{det}{Proportion of recurrent points forming diagonal + \item{DET}{Proportion of recurrent points forming diagonal line structures.} - \item{nrline}{The total number of lines in the recurrent plot} - \item{maxline}{The length of the longest diagonal line + \item{NRLINE}{The total number of lines in the recurrent plot} + \item{maxL}{The length of the longest diagonal line segment in the plot, excluding the main diagonal} - \item{meanline}{The average length of line structures} - \item{entropy}{Shannon information entropy of + \item{L}{The average length of line structures} + \item{ENTR}{Shannon information entropy of diagonal line lengths longer than the minimum length} - \item{relEntropy}{Entropy measure normalized by the number of + \item{rENTR}{Entropy measure normalized by the number of lines observed in the plot. Handy to compare across contexts and conditions} - \item{lam}{Proportion of recurrent points forming vertical + \item{LAM}{Proportion of recurrent points forming vertical line structures} - \item{tt}{The average length of vertical line structures} + \item{TT}{The average length of vertical line structures} } \author{Moreno I. Coco (moreno.cocoi@gmail.com)} @@ -83,10 +92,10 @@ mindiagline, minvertline, tw, whiteline, recpt) version provided by Rick Dale, and created during the Non-Linear Methods for Psychological Science summer school held at the University of Cincinnati - in 2011} + in 2012} -\seealso{ \code{\link{tt}}, \code{\link{spdiags}}, -\code{\link{simts}}} +\seealso{ \code{\link{tt}}, \code{\link{checkts}}, + \code{\link{spdiags}}, \code{\link{simts}}} \examples{ @@ -98,12 +107,15 @@ ts1 = tS[1,]; ts2 = tS[2,] ## (e.g., RDts1, RDts2) ## initialize the parameters -delay = 1; embed = 1 ; rescale = 1; radius = 0.00001; -normalize = 0; minvertline = 2; mindiagline = 2; -tw = 0; whiteline = FALSE; recpt = FALSE -ans = crqa(ts1, ts2, delay, embed, rescale, radius, -normalize, minvertline, mindiagline, tw, whiteline,recpt) +delay = 1; embed = 1; rescale = 1; radius = 0.001; +normalize = 0; mindiagline = 2; minvertline = 2; +tw = 0; whiteline = FALSE; recpt = FALSE; side = "both" +checkl = list(do = FALSE, thrshd = 3, datatype = "categorical", + pad = TRUE) + +ans = crqa(ts2, ts1, delay, embed, rescale, radius, normalize, +mindiagline, minvertline, tw, whiteline, recpt, side, checkl) print(ans[1:9]) ## last argument of list is the cross-recurrence plot RP = ans$RP ## take out RP diff --git a/man/optimizeParam.Rd b/man/optimizeParam.Rd index fc6e0e5..11a11e7 100644 --- a/man/optimizeParam.Rd +++ b/man/optimizeParam.Rd @@ -14,7 +14,7 @@ } \usage{ -optimizeParam(ts1, ts2, par) +optimizeParam(ts1, ts2, par, min.rec, max.rec) } \arguments{ @@ -27,17 +27,20 @@ optimizeParam(ts1, ts2, par) the two series. steps = a sequence of points (e.g., seq(1, 10, 1)) used to look ahead local minima. - cut.del = a sequence of points referring to - the delays evaluated when mutual information - between the two serie is estimated. + fnnpercent = the percentage of information gained, + relative to the first embedding dimension, when higher embeddings + are considered. radiusspan = a constant setting the granularity of radius unit to explore, relative to the standard deviation of the distance between the two series. (Larger value = smaller units) radiussample = the number of equally spaced units of the radius to explore. (Larger value = more units = - computationally more expensive) - - } + computationally more expensive)} + \item{min.rec}{The minimum value of recurrence accepted. + Default = 2} + \item{max.rec}{The maximum value of recurrence accepted. + Default = 5} + } \details{ @@ -73,12 +76,13 @@ optimizeParam(ts1, ts2, par) } \note{As \code{optimizeParam} uses \code{crqa} to estimate - the parameters: the additional arguments \code{normalize, + the parameters: the additional arguments normalize, rescale, mindiagline, minvertline, whiteline, - recpt} should be supplied in the par list. + recpt should be supplied in the par list. Set up relatively large radiusspan (e.g. 100), and relatively small radiussample (e.g., 10), - for a decent coverage of radius values.} + for a decent coverage of radius values. + } \value{ @@ -97,7 +101,8 @@ optimizeParam(ts1, ts2, par) } -\author{Moreno I. Coco (moreno.cocoi@gmail.com)} +\author{Moreno I. Coco (moreno.cocoi@gmail.com) and James A. Dixon +(james.dixon@uconn.edu)} \seealso{ \code{\link{crqa}}, \code{\link{wincrqa}}} @@ -105,18 +110,19 @@ optimizeParam(ts1, ts2, par) ## initialize the parameters -par = list(lgM = 20, steps = seq(1, 6, 1), -cut.del = seq(1, 40,1), +par = list(lgM = 20, steps = seq(1, 6, 1), radiusspan = 100, radiussample = 40, normalize = 0, rescale = 1, mindiagline = 2, minvertline = 2, tw = 0, whiteline = FALSE, -recpt = FALSE) +recpt = FALSE, fnnpercent = 10) + ## generate two random uniform series + ts1 = runif(100) ts2 = runif(100) -ans = optimizeParam(ts1,ts2,par) +ans = optimizeParam(ts1,ts2,par, min.rec = 2, max.rec = 5) print(ans) ## utilize leftmov, rightmov for an application to real-data diff --git a/man/runcrqa.Rd b/man/runcrqa.Rd index e7c99c8..19ea9a9 100644 --- a/man/runcrqa.Rd +++ b/man/runcrqa.Rd @@ -172,7 +172,8 @@ data(crqa) ###Cross-recurrence diagonal profile par = list(type = 1, ws = 100, method = "profile", - datatype = "continuous", thrshd = 8, radius = 2) + datatype = "categorical", thrshd = 8, radius = 2, + pad = FALSE) ans = runcrqa(RDts1, RDts2, par) @@ -181,8 +182,9 @@ profile = ans$profile; maxrec = ans$maxrec; maxlag = ans$maxlag #################################################### ###Windowed cross-recurrence profile -par = list(type = 1, step = 20, windowsize = 100, lagwidth = 40, - method = "window", datatype = "categorical", thrshd = 8) +par = list(type = 1, windowstep = 20, windowsize = 100, + method = "window", datatype = "categorical", thrshd = 8, + pad = FALSE) ans = runcrqa(RDts1, RDts2, par) @@ -199,4 +201,4 @@ res = runcrqa(RDts1, RDts2, par) res[1:9] -} \ No newline at end of file +} diff --git a/man/tt.Rd b/man/tt.Rd index acc237b..cc9c10a 100644 --- a/man/tt.Rd +++ b/man/tt.Rd @@ -38,7 +38,7 @@ \author{Moreno I. Coco (moreno.cocoi@gmail.com)} -\details{This function is based on original MATLAB code +\details{This function is based on MATLAB code written by Norbert Marwan, available in \code{crptoolbox}} \examples{ diff --git a/man/wincrqa.Rd b/man/wincrqa.Rd index 4da9d1a..1913916 100644 --- a/man/wincrqa.Rd +++ b/man/wincrqa.Rd @@ -14,19 +14,18 @@ \usage{ -wincrqa(x, y, step, windowsize, lagwidth, delay, embed, rescale, -radius, normalize, mindiagline, minvertline, tw, whiteline, recpt) +wincrqa(ts1, ts2, windowstep, windowsize, delay, embed, rescale, +radius, normalize, mindiagline, minvertline, tw, whiteline, trend) } \arguments{ - \item{x}{First time-series} - \item{y}{Second time-series} - \item{step}{Interval by which the window is moved.} + \item{ts1}{First time-series} + \item{ts2}{Second time-series} + \item{windowstep}{Interval by which the window is moved.} \item{windowsize}{The size of the window} - \item{lagwidth}{The number of delays to be considered} \item{delay}{The delay unit by which the series are lagged.} \item{embed}{The number of embedding dimension for phase-reconstruction, i.e., the lag intervals.} @@ -47,9 +46,8 @@ radius, normalize, mindiagline, minvertline, tw, whiteline, recpt) \item{tw}{The size of the Theiler window} \item{whiteline}{A logical flag to calculate (TRUE) or not (FALSE) empty vertical lines.} - \item{recpt}{A logical flag indicating whether - measures of cross-recurrence are calculated directly - from a recurrent plot (TRUE) or not (FALSE).} + \item{trend}{A logical flag indicating whether + the TREND should be computed} } @@ -69,7 +67,14 @@ radius, normalize, mindiagline, minvertline, tw, whiteline, recpt) \note{If no-recurrence is found in a window, that window will not be saved, and a message - about it will be warned.} + about it will be warned. + TREND is implemented following a solution proposed + by Norbert Marwan, and translated here in R, + for those who have asked him. + He, however warns that this measure might strongly depend + on the chosen settings to calculate crq. + Relying on such measure can, therefore, produce misleading results. + } \seealso{\code{\link{crqa}}} @@ -77,24 +82,25 @@ radius, normalize, mindiagline, minvertline, tw, whiteline, recpt) \examples{ ## simulate two dichotomous series -tS = simts(0.25, 0.05, 0.2, 0.2, 0.25, 200) +tS = simts(0.25, 0.05, 0.2, 0.2, 0.25, 100) ts1 = tS[1,]; ts2 = tS[2,] ## check data(crqa) for alternative data ## (e.g., RDts1, RDts2) -step = 10; windowsize = 100; lagwidth = 20; +windowstep = 10; windowsize = 50; delay = 1; embed = 1 ; rescale = 1; radius = 0.00001; normalize = 0; minvertline = 2; mindiagline = 2; -tw = 0; whiteline = FALSE; recpt = FALSE +tw = 0; whiteline = FALSE; trend = TRUE; -## it returns matrix with all measures for the -## different windows where values are found +## it returns a list with: +## [[1]] the measures for the different windows where values are found +## [[2]] the trend over time. -res = wincrqa(ts1, ts2, step, windowsize, -lagwidth, delay, embed, rescale, radius, normalize, -mindiagline, minvertline, tw, whiteline, recpt) +res = wincrqa(ts1, ts2, windowstep, windowsize, +delay, embed, rescale, radius, normalize, mindiagline, +minvertline, tw, whiteline, trend) str(res)