Danny Kaplan’s nonlinear time series functions

The Mathematica functions below are translations of D. Kaplan’s Matlab code available at

Off[General::spell]

Acknowledgments

Firstly, I obviously have to thank Daniel Kaplan for allowing me to publish the translation of his code, and for helping me when I got stuck. Mathew Dafilis of Swinburne University in Melbourne also helped me fixed some remaining bugs.

Packages to Load in

Needs["Statistics`MultiDescriptiveStatistics`"]

Required for plotACF

Needs["Statistics`DescriptiveStatistics`"]

Required for theilerQ, crinkleStatistic, nonlinearPredictionError

Needs["Statistics`DataManipulation`"]

Required for the various box* functions and compHistogram

Needs["Utilities`FilterOptions`"]

Required for kaplanScatterPlot, plotDE  and plotACF

Needs["Graphics`Graphics`"]

Required for compHistogram

Needs["Graphics`Graphics3D`"]

Required for kaplanScatterPlot

Needs["LinearAlgebra`MatrixManipulation`"]

Required for midCast,timeReversalAssymetryStatistic

Needs["LinearAlgebra`Orthogonalization`"]

Needs["Statistics`ContinuousDistributions`"]

Required for model1, model1realization,  model2, and model2realization.

$Packages

Functions not implemented

A number of the functions implemented in Kaplan's set of M-files are already implemented in the standard distribution of Mathematica, or can be built up from basic Mathematica functions.  These include Kaplan's  gausspdf.m (see "Statistics`ContinuousDistributions`"), rho.m (See "Statistics`MultiDescriptiveStatistics`"), globelin.m (see "Statistics`LinearRegression`"), bootstr.m and comphist.m. Similarly, repeatp.m is probably easiest to achieve using a function like NestList, FoldList or Table, and Mapping the desired function onto the resultant list, all of which can be done on a single line. scatter3movie.m and scattermovie.m  are best implemented using Mathematica's front end commands for generating movies in various formats.
surrstat.m just maps a particular statistic onto surrogate data. Rather than write a function to do this specifically, it seems that it would be best to use nameOfFunction[fftSurrogate[data]]  directly.
fditer.m just iterates on a function and returns a list of the iterates.  May as well use NestList[].
squid1.m and squid2.m are only half-implemented because they seem to refer to external code, which does not appear to be included in the download. I didn't do squidper.m for the same reason.
nearneib.m doesn't appear to do anything.
xfershow.m relies on function that are not defined in other code in the package, namely coherence, that I don't know how to implement.
polyshow.m can be replicated using Mathematica's built-in Fit[] function and the plots of your choice. 
Similarly, threeplot.m can be implemented using built-in Plot functions (check out the documentation on Statistics`DataManipulation`BinCounts and BarChart to get histograms with a minimum of fuss.


 I have also decided not to implement the Lorenz ODE, since this is in the Mathematica Book, section 3.9.7, (this link might only work on Macs, as it is path dependent).
NDSolve[ {x'[t] == -3 (x[t] - y[t]),  y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],   z'[t] == x[t] y[t] - z[t],   x[0] == z[0] == 0, y[0] == 1},   {x, y, z}, {t, 0, 20}, MaxSteps->3000 ]

Code by Luci Ellis

simpleLinRegress[ ]: A Simple Linear Regression function built for speed not features

simpleLinRegress::usage = 
      "simpleLinRegress[data_MatrixQ] takes a rectangular matrix of real \
numbers and performs a linear regression or weighted linear regression of the \
first column on the other columns. The output is a non-rectangular list of \
three lists -- the estimated coefficients, the fitted values and the \
residuals."

GLSWeights::usage  =  
      "GLSWeights is an option to simpleLinRegress that takes either the \
value Automatic, corresponding to OLS, or a list of real numbers \
corresponding to the weights to be used in a WLS procedure."

IncludeConstant::usage  =  
      "IncludeConstant is an option to simpleLinRegress that indicates \
whether or not a constant should be added to the data matrix."

simpleLinRegress::wtnum = 
      "Either the vector of weights is not a vector of real numbers with \
the same number of element as the data sample size, or you have tried to \
apply WLS to a model with a single explanator and no constant."

simpleLinRegress::notnum  =  
      "Some of the data supplied is non-numeric."

Options[simpleLinRegress] = {GLSWeights -> Automatic, 
        IncludeConstant -> True}

SetOptions[simpleLinRegress, GLSWeights -> Automatic, 
      IncludeConstant -> True]

\(simpleLinRegress[data_?MatrixQ, \ opts___?OptionQ] := \ 
      If[And\ @@\ \((NumericQ\ /@\ Flatten[data])\), 
		
        Module[{𝓌, 𝒸}, 
          𝓌\  = \ 
            \(GLSWeights\  /. {opts}\)\  /. Options[simpleLinRegress]; 
		
          𝒸\  = \ 
            \(IncludeConstant\  /. {opts}\)\  /. Options[simpleLinRegress]; 

          		With[{
              𝓎 = 
                If[Not[𝒸] && \ 
                    \(Dimensions[data]\)[[2
                        ]] == 2, Flatten[First/@data], 
                  First\ /@\ data], 
							\ \ 
              𝓍\  = \ 
                If[𝒸, 
                  Transpose[
                    Join[{Table[1, {Length[data]}]}, 
                      Transpose[Rest\ /@\ data]]], 
                  If[\(Dimensions[Rest\ /@data]\)[[2
                        ]] == 1, Flatten[Rest\ /@\ data], 
                    Transpose[Rest\ /@\ data]]], 𝓃 = Length[data]}, 
            
			Module[{betas, yhat, \ resids}, 
				
              Which[𝓌 === Automatic, 
                \(betas\  = \ 
                  If[Not[𝒸] && 
                      \(Dimensions[Rest\ /@data]\)[[2
                          ]] == 1, {
                      \((1/\((𝓍.𝓍)\))\)*
                        \((𝓍.𝓎)\)}, 
                    Flatten[
                      Inverse[\((Transpose[\ 𝓍].\ 𝓍)\)].
                        \((Transpose[\ 𝓍].𝓎)\)]]; \), 

                					
                \((Length[𝓌] == 𝓃\ )\) && 
                  \((Length[𝓌] > 1\ )\) && \ 
                  \((And\ @@\ \((NumericQ\ /@\ 𝓌)\))\), \ 
                \(betas\  = 
                  Inverse[
                      Transpose[\((𝓌*𝓍)\)].
                        \((𝓍)\)].
                    \((\((Transpose[𝓌*𝓍])\).
                        \((𝓎)\))\); \), True, \ 
                Message[simpleLinRegress::wtnum]; \ $Failed]; 
				
              yhat\  = 
                If[Not[𝒸] && \ 
                    \(Dimensions[data]\)[[2
                        ]] == 2, 
                  𝓍*betas[[1]], 
                  \ 𝓍.betas]; 
				
              resids\  = \ 𝓎 - yhat; 

              				{betas, yhat, resids}]]], 
        Message[simpleLinRegress::notnum]]\)

mHistory

mHistory::usage = 
      "mHistory[data,m] constructs a m-embedding of the data series data \
with lag 1. mHistory[data,m,l] constructs an m-embedding with lag l. This \
function replicates the functionality of a Gauss procedure lagembed, written \
by Daniel Kaplan, but was written by Luci Ellis prior to her seeing his code.\
"

mHistory[data_?VectorQ, 1] :=  data

mHistory[data_?VectorQ, m_Integer] :=  
      Drop[Transpose[RotateLeft[data, #] & /@ Range[0, m - 1]], 1 - m]

\(mHistory[data_?VectorQ, m_Integer, lag_Integer] /; 
        Length[data] > \((m + lag + 1)\) := 
	\ 
      \(Drop[Transpose[\(RotateLeft[data, #]\ &\)\ /@\ Range[0, m - 1]], 
          1 - m - lag]\)[[Range[1, Length[data] - m - lag, 
          lag]]]\)

mHistory[data_?MatrixQ, 1] :=  data

mHistory[data_?MatrixQ, m_Integer] :=  
      Drop[Transpose[RotateLeft[data, #] & /@ Range[0, m - 1]], 1 - m]

\(mHistory[data_?MatrixQ, m_Integer, lag_Integer] /; 
        Length[data] > \((m + lag + 1)\) := 
	\ 
      \(Drop[Transpose[\(RotateLeft[data, #]\ &\)\ /@\ Range[0, m - 1]], 
          1 - m - lag]\)[[Range[1, Length[data] - m - lag, 
          lag]]]\)

Code by Danny Kaplan translated from Matlab to Mathematica by Luci Ellis

plotACF

plotACF::usage="plotACF plots the autocovariance function up to a given lag, showing the 95 percent confidence intervals for white noise.
This version is based on Matlab code by Daniel Kaplan."

\(plotACF[data : {__?NumericQ}, maxlag_Integer, opts___Rule] /; 
        \((maxlag < \((Length[data]/2)\))\) := 

      With[{n = Length[data]}, 
        With[{confint = Table[2/Sqrt[n], {maxlag}]}, 
			
          Show[ListPlot[
              \(Correlation[Drop[data, #], Drop[data, \(-#\)]]&\)\ /@\ 
                Range[maxlag], DisplayFunction -> Identity], 
			
            ListPlot[confint, PlotJoined -> True, 
              PlotStyle -> Dashing[{0.02, \ 0.02}], 
              DisplayFunction -> Identity], 
			
            ListPlot[\(-confint\), PlotJoined -> True, 
              PlotStyle -> Dashing[{0.02, \ 0.02}], 
              DisplayFunction -> Identity], 
            FilterOptions[ListPlot, opts], 
            DisplayFunction -> $DisplayFunction]\ ]]\)

abcalc

abcalc::usage="abcalc[alpha0, beta0, alpha1, beta1, p0, p1] calculates a and b for feedback control of one-dimensional maps for the definitions of the arguments, see the Squid Lab instructions. This version is based on Matlab code by Daniel Kaplan."

\(abcalc[α0_, \ β0_, \ α1_, \ β1_, \ p0_, \ p1_]
         := 
	With[{
          γ\  = \ \((β1 - β0)\)/\((p1 - p0)\)}, 
        
{\ \(-α0\)\ /\((γ\  - α0\ *\ γ)\), \ 
          \((α0*β0\  + γ*p0\  - \ α0*γ*p0)
              \)/\((γ\  - α0*γ)\)}]\)

ampsurr

ampsurr::usage = "ampsurr[x,y,seed] generates amplitude adjusted surrogate data. x is the data, y is data from a target histogram, or an Integer.
n == 1  ==> N(0,1)
n == 2  ==> Uniform(0,1) --- this is deterministic
n == 3  ==> Triangular(0,2)
This version is based on Matlab code by Daniel Kaplan."

\(ampsurr[x_?VectorQ, y_:  1, seed_] := 

      Module[{}, SeedRandom[seed]; 
\ 
        histxfrm[fftSurrogate[\ histxfrm[x, y, seed], seed], x]]\)

\(ampsurr[x_?VectorQ, y_:  1] := 
\ 
      histxfrm[fftSurrogate[\ histxfrm[x, y]], x]\)

ARResids

ARResids::usage="ARResids[data, porder]
calculates the residuals of data from the specified AR model. This function is based on Matlab code by Daniel Kaplan and uses the simpleLinRegress quick-and-dirty OLS function by Luci Ellis.  The difference between this version and the Kaplan version is that this one just takes the order (p) of the AR model to be estimated, estimates it by OLS and returns the residuals."

\(ARResids[data : {__?NumericQ}, p_Integer] /; 
        \((p\  < Length[data])\) := 
	
      With[{olsdata = 
            Transpose[
              \(Drop[#, p]&\)\ /@
                \((\(RotateRight[data, #]&\)\ /@\ Range[0, p])\)]}, 
        Last[simpleLinRegress[olsdata]]]\)

bincountsInformation

bincountsInformation::usage="bincountsInformation[bincounts] calculates the information given a probability distribution in the form of bincounts.  
Example usage:  (s,t) = hist(data); inform(s) r is normalized to the information in the uniform distribution with length(bincnts).
This version is based on Matlab code by Daniel Kaplan."

\(bincountsInformation[bincnts_List] := 

      With[{counts = \ Select[\ bincnts\ , # > \ 0\ &]}, 
		
        With[{normcounts = counts/\((Plus\ @@\ counts)\)}, 

          \(-\((\ Plus\ @@\((normcounts*Log[normcounts])\))\)\)/
            Log[Length[bincnts]]\ ]]\)

bootstrap

bootstrap::usage="bootstrap[data,stat,n] estimates the standard error of a statistic using bootstrapping.
data is a matrix containing the data --- sampling is done row-by-row
stat is the statistic to use. 
n is the number of bootstrap samples to take.
This version is based on Matlab code by Daniel Kaplan."

\(bootstrap[data, \ stat_Symbol, \ n_Integer] := 
	
      With[{s = Dimensions[data]}, 
        With[{len = Max[s], 
            datatouse = 
              If[s[[1]] < 
                  s[[2]], 
                Transpose[data], data]}, 
	
          Table[Evaluate[
              stat\ @@\ 
                datatouse
                  [[Ceiling[len*Table[Random[], {len}]]\ 
                  ]]], {n}]\ ]]\)

boxCount and boxDim

boxCount::usage="boxCount[ data, nbins] counts the number of boxes that a matrix occupies. Uses the standard package Statistics`DataManipution.
This version is based on Matlab code by Daniel Kaplan."

\(boxCount[data_?MatrixQ, nbins_Integer] := 
	
      With[{mx = Ceiling[Max[data]], mn = Floor[Min[data]]}, 
        With[{bincnts = BinCounts[data, {mn, mx, \((mx - mn)\)/nbins}]}, 
          Length[Select[bincnts, #\  != \ 0\ &]]/nbins]]\)

An equivalent way to get the same answer.  It remains to be tested which is more efficient.

\(boxCount2[data_?MatrixQ, nbins_Integer] := 
	
      With[{mx = Ceiling[Max[data]], mn = Floor[Min[data]]}, 
        With[{bincnts = 
              BinCounts[Flatten[data], {mn, mx, \((mx - mn)\)/nbins}]}, 
          \((nbins - Length[Cases[bincnts, 0]])\)/nbins]]\)

boxDim::usage="boxDim[data] returns two columns, 
 log of the box resolution and log of the number of occupied boxes at that resolution.
This version is based on Matlab code by Daniel Kaplan."

\(boxDim[data_?MatrixQ, n_Integer] := 
	
      Table[{Log[i], Log[i*\ boxCount[data] + 0.0000001]}, {i, n}]\)

compHistogram (comphist.m)

compHistogram::usage = 
      "compHistogram[data, tname, numbins] overplots a histogram and a \
normal probability
 distribution. This version is based on Matlab code by Daniel Kaplan. The \
x-axis labels are only displayed to one decimal place."

\(compHistogram[data_?VectorQ, \ tname_String, numbins_Integer, 
        opts___Rule]\  := 
	
      With[{μ = Mean[data], σ = StandardDeviation[data], 
          α = Min[data], β = Max[data], 
          𝓃 = Length[data]}, 
		
        With[{j = Abs[\((\((β - α)\)/\((numbins)\))\)]}, 
		
          Show[{\ BarChart[
                Transpose[{
                    BinCounts[data, 
                      \ {α, \ β, \ 
                        \((\((β - α)\)/\((numbins)\))\)}]\ , 
                    Round[10\ 
                          Rest[Range[α, β, \ j] - \((j/2)\)\ ]]/
                      10.}], DisplayFunction -> Identity], \ 

              ListPlot[
                Table[𝓃\ *
                      \((CDF[NormalDistribution[μ, σ], x + j])\) - 
                    𝓃\ *
                      CDF[NormalDistribution[μ, σ], x], {x, 
                    α, α + j\ numbins, j\ }], 
                PlotRange -> All, PlotJoined -> True, 
                DisplayFunction -> Identity]}, PlotLabel -> tname, 
            DisplayFunction -> $DisplayFunction, 
            FilterOptions[ListPlot, opts]]\ ]]\)

evalMap

evalMap::usage="evalMap[zz,targs,k] finds the values of a locally linear map and a locally constant map for targs = f(zz)
zz - embedding of data
targs -- the images of zz under the dynamics --- the map will be fit to this.
k - number of neigbors to use in constructing the map
RETURNED VALUES
lin -- values of the locally linear map at each point in zz
const -- values of the locally contant map at each point. This version is based on Matlab code by Daniel Kaplan."

evalMap[zz_,targs_,k_]:=
	With[{len =Length[zz]},
Transpose[Table[localLinearFitMap[zz, zz[[j]], k, targs, j]
	(* exclude the current point j from the fitting *),{j,len}]][[{3,4}]] ]

filterSpecification (filtspec.m)

gaussTruth::usage = 
      "gaussTruth is just a one-liner that replicates the functionality of \
the way Gauss and Matlab do greater than / less than tests on lists, \
returning (0,1) instead of True/False the way Mathematica does. However, note \
that if the value tested is exactly equal to 0, gaussTruth returns 1/2 so you \
might have to explicitly wrap a Sign[] function around the gaussTruth \
command. As an example, to replicate vector.< x , use \
gaussTruth[x-vector]."

gaussTruth[x_] :=  1/2  +  Sign[x]/2

filterSpecification::usage = 
      "filterSpecification[ cut, len, width] is a 
function for specifying a low-pass filter in the fourier domain
cut - cutoff frequency in terms of the sampling frequency.
len - length of the data to be filtered.
width - width of the transition zone.
This version is based on Matlab code by Daniel Kaplan."

\(filterSpecification[\ cut_, \ len_, \ width_] /; 

        \((\((	cut\  <= \ 0.5)\)\  && \ \((cut\  + \ width\  <= \ 0.5)\))
          \)\) :=
\(Module[{t, τ}, 	
        With[{high1 = len/2, 
            high2 = If[Mod[len, 2] == 0, len/2 - 1, \ Floor[len/2]\ ]}, 
		
          With[{q\  = \ {\ 0, high1\ high2, \(-1\), 1\ }/len}, 
			
            τ = \ \((cut\  + \ width\  - q)\)/\ width; 
			
            t = If[width == 0, 0*q, 
                τ\ *gaussTruth[0.9999 - τ]*
                  gaussTruth[τ - 0.0001]]; 

            t\  + Sign[gaussTruth[cut - q]]\ \ ]]]\)

findNeighbour (findneib.m)

findNeighbour::usage="findNeighbour[data_, pt, k, r ] finds the nearest neighbors to pt in z
data -- matrix of points, 1 per row
pt -- vector of a single point 
k -- number to find
r (optional )-- if specified, find all neighbors closer than this
inds -- indices of the closest points to pt
dist -- corresponding distances from pt
 This version is based on Matlab code by Daniel Kaplan."

\(findNeighbour[data_?MatrixQ, \ pt_?VectorQ, \ k_, \ r_:  0\ ] /; 
        \((\(Dimensions[data]\)[[2]]\ 
             == \ Length[pt])\)\  := 

      With[{ds\  = \ 
            Sort[Transpose[{oneDist[data, pt], Range[Length[data]]}\ ]]\ }, 

        Reverse\ /@\ 
          \((If[r == 0, Transpose[Take[ds, k]], 
              Transpose[
                Select[ds, #[[1]]\  < \ r
                      \ &]]\ ])\)]\)

\(findNeighbour[data_?VectorQ, \ pt_?NumericQ, \ k_, \ r_:  0\ ] := 

      With[{ds\  = \ 
            Sort[Transpose[{oneDist[data, pt], Range[Length[data]]}\ ]]\ }, 

        Reverse\ /@\ 
          \((If[r == 0, Transpose[Take[ds, k]], 
              Transpose[
                Select[ds, #[[1]]\  < \ r
                      \ &]]\ ])\)]\)

fftFilter

fftFilter::usage="fftFilter filters in the spectral domain. The filter is usually constructed using filtSpec.
This version is based on Matlab code by Daniel Kaplan."

\(fftFilter[data : {__?NumericQ}, filter_?VectorQ] /; 
        \((Length[data] == Length[filter]\ )\) := 

      Re[InverseFourier[Fourier[data]*filter]]\)

fftSurrogate

fftSurrogate::usage="fftSurrogate calculates phase randomized surrogate data from a time series x, with an optional specified seed. This version is based on Matlab code by Daniel Kaplan, with bugfixes by Mathew Dafilis."

fftSurrogate[x_?VectorQ, seed_?NumericQ]:=
Module[{ph, iph},	With[{z = Fourier[x]},
		SeedRandom[seed];
(* randomize the phases and replace 2 Pi with some
                         other phase since 2 Pi  ==  0 *)
				ph  =  
                  (2*Pi*
                        Table[Random[], {1, Ceiling[Length[z]/2  -  1]}] )
                       //.  N[2 Pi]  -> Random[Real, {0, 2  Pi}];
	\(iph\  = 
                  Flatten[
                    If[\ Mod[Length[z], 2]\  == \ 0\ , 

                      \ {\ 0, \ ph, \ 0, \ \(-Reverse[ph]\)}, 
                      \ {\ 0, \ ph, \ \ \(-Reverse[ph]\)\ }]]\);
Re[InverseFourier[ z* Exp[I*iph]]] ]]

fftSurrogate[x_?VectorQ]:=
Module[{ph, iph},	With[{z = Fourier[x]},
  
				\(ph\  = \ 
                \((2*Pi*Table[Random[], {1, Ceiling[Length[z]/2\  - \ 1]}]\ )
                    \) //. \ N[2\ Pi]\  -> Random[Real, {0, 2  Pi}]; 
	
              iph\  = Flatten[
                  If[\ Mod[Length[z], 2]\  == \ 0\ , 

                    \ {\ 0, \ ph, \ 0, \ \(-Reverse[ph]\)}, 
                    \ {\ 0, \ ph, \ \(-Reverse[ph]\)\ }]]; 

              Re[InverseFourier[\ z*\ Exp[I*iph]]]\) ]]

histogramTransformation (histxform.m)

?? histosort

histosort::usage = 
      "histosort[x,y] is an auxiliary function that is called by \
histogramTransformation, to sort the data in y into the same order implied by \
the data in x. Both x and y should be N*2 matrices with the data in the first \
column and integers in the second column."

\(histosort[xs : {{_, _Integer}..}, ys : {{_, _Integer}..}] /; 
        \((\(Dimensions[xs]\)[[1]] == 
            \(Dimensions[ys]\)[[1]])
          \) := 
	
      With[{originalindices = 
            Flatten[Rest\ /@\ 
                \((Sort[Transpose[{\((Rest\ /@\ xs)\), Range[Length[xs]]\ }]])
                  \)]}, 
		
        First\ /@\ 
          ys[[originalindices]]\ ]\)

Note that, strictly speaking, I could have just piled all the calculations to obtain the indices inside the part specification: 
ys[[Sort[Transpose[{(Rest /@ xs),Range[Length[xs]] }]]]], without using a local constant and a With statement.
I decided to do it the way I did for ease of understanding.

histogramTransformation::usage = 
      "histogramTransformation[x,y] transforms x so that it has the same \
histogram as y while maintaining the order of the ranks in x. x and y should \
be single column vectors.
histogramTransformation[x,dist,seed] transforms x to have a histogram of the \
type specified by the distribution dist, with (optional) seed s.
This version is based on Matlab code by Daniel Kaplan."

\(histogramTransformation[x_?VectorQ, y_?VectorQ] /; 
        \((Length[x] > 0\  && Length[y] > 0)\) := 
	
      With[{xrange = Range[Length[x]]}, 
		
        With[{xs = Sort[Transpose[{x, xrange}]], 
            ys = Sort[Transpose[{y, xrange}]]}, histosort[xs, ys]\ ]]\)

\(histogramTransformation[x_?VectorQ, dist_] /; \((Length[x] > 0)\) := 

      	With[{y = RandomArray[dist, Length[x]], xrange = Range[Length[x]]}, 

        		With[{xs = Sort[Transpose[{x, xrange}]], 
            ys = Sort[Transpose[{y, xrange}]]}, histosort[xs, ys]]]\)

\(histogramTransformation[x_?VectorQ, dist_, seed_?NumericQ] /; 
        \((Length[x] > 0)\) := 
	
      Module[{}, SeedRandom[seed], histogramTransformation[x, dist]]\)

localMap

localMap::usage="localMap[z, pt, k, exclude] constructs the
indices of neighboring points from z, a matrix of embedded data; pt, the point at which to find the neighbors; and k, the number of nearest neighbors to use exclude (optional) exclude this point from the preimages.  This version is based on Matlab code by Daniel Kaplan."

\(localMap[z_?MatrixQ, \ pt_?VectorQ, \ k_, \ exclude_Integer] := 

      Transpose[Delete[\ Transpose[\ findNeighbour[z, pt, k + 1]], exclude]]
        \)

\(localMap[z_?VectorQ, \ pt_?NumericQ, \ k_, \ exclude_Integer] := 

      Transpose[Delete[\ Transpose[\ findNeighbour[z, pt, k + 1]], exclude]]
        \)

localMap[z_?MatrixQ,  pt_?VectorQ,  k_] := findNeighbour[z, pt, k]

localMap[z_?VectorQ,  pt_?NumericQ,  k_] := 
      findNeighbour[z, pt, k]

localLinearFitMap (locallin.m) (This doesn't quite work yet -- The LinearSolve seems to come up with no solution!)

The reason is that Transpose[Join[{Table[1.,{Length[f]}]},Transpose[ p]]] (p with 1's as the first column) has the same length as f, but it's not square.  That is, the system is underdetermined (x equations in y unknowns, where x > y by a long way!).  I'm not sure why this happens yet, since the function depends on so many other functions it's hard to trace through.

localLinearFitMap::usage="localLinearFitMap[z,pt,k,targs,exclude] finds a locally linear map for z(t+1) based on embedded values where z is embedded data; pt is the point at which to construct the map; k is how many neighbors to use; targs -- the actual target values in the time series; exclude is a point to be excluded from the map. The output is {a0,a1,val,mn} where a0 is the constant in map; a1 is the vector by which to multiply pt; val is the value of the map evaluated at pt and mn is value of a locally constant map. This version is based on Matlab code by Daniel Kaplan."

localLinearFitMap[z_, pt_, k_, targs_, exclude_]:=
	Module[{f, p, foo},
With[{lm = (localMap[z, pt, k, exclude])},
			With[{inds = lm[[2]]},
					blah = inds;
f =  (
                        targs[[inds]])
                        ;
p  =  
                      z[[inds]];
  
foo  = 
                      LinearSolve[
                        Transpose[
                          Join[{Table[1., {Length[f]}]}, Transpose[ p]]], f]
                        ;
{ First[foo],  Rest[foo],   
                      First[foo]  +  pt*Rest[foo],  Mean[f]}]]]

localLinearFitMap[z_, pt_, k_, targs_]:=
		Module[{f, p, foo},With[{lm = (localMap[z, pt, k])},With[{inds = lm[[2]]},
								blah = inds;
f =   
                      (targs[[inds]])
                        ;
p  =  
                      z[[inds]];
  \(

                    
                    foo\  = 
                      LinearSolve[
                        Transpose[
                          Join[{Table[1., {Length[f]}]}, Transpose[\ p]]], 
                        f]; 
{\ First[foo], \ Rest[foo], \ 
                      First[foo]\  + \ pt*\ Rest[foo], \ Mean[f]}\)]]]

mapControl

mapControl::usage="mapControl[x0,map,a,b,n,mu0,xon,xsize] iterates a map using linear feedback control.
x0 -- initial condition
map -- the map to iterate map[x,param]
a, b -- parameters in the feedback control
n -- number of iterates
μ0 -- default value of parameter
xon, xsize -- turn on control when x is in xon+-xsize.
This version is based on Matlab code by Daniel Kaplan. It is similar to a standard NestList except that the second parameter passed to the map varies according to the previous value of of iterated function."

\(mapControl[x0_, thismap_Symbol, a_, b_, n_Integer, μ0_, 
        xon_Integer, xsize_Integer] := 

      Rest[NestList[\ 
          \((thismap
                \ @@{#, 
                  If[Abs[# - xon]\  < \ xsize, \ a*# + \ b, \ \ μ0]})
              \)\ &, x0, n]]\)

Applying a Symbol that is the name of a map works like this:

f @@ {a, b}

midCast

midCast::usage = 
      "midCast[z] is similar to getImage[] but for midcasting.
This version is based on Matlab code by Daniel Kaplan. It requires the \
MatrixManipulation package from the LinearAlgebra group."

\(midCast[z_?MatrixQ] /; 
        \((\(Dimensions[z]\)[[2]] > 2)
          \) := 
With[{
          col\  = \ 
            Round[.51\  + \ 
                \(Dimensions[z]\)[[2]]/
                  2]}, 
{\ \(Delete[#, col]&\)\ z\ , \ 
          TakeColumns[z, {col, col}]}\ ]\)

model1

model1::usage = 
      "model1[α,β,ν,x0] runs model 1 dynamics and \
measurements given the  noise inputs ν. x0 is the initial condition. This \
version is based on Matlab code by Daniel Kaplan."

\(model1[α_, β_, ν_?VectorQ, x0_] := 
	
      ν\  + 
        Rest[NestList[α\ *\ #\  + \ β\ &, x0, Length[ν]]]\)

\(model1[α_, β_, ν_?VectorQ] := 
	
      With[{x0\  = \ β/\((1 - α)\)}, 
        model1[α, β, ν, x0]\ ]\)

\(model1[α_, β_] := 
	
      With[{x0\  = \ β/\((1 - α)\), 
          ν = \ RandomArray[NormalDistribution[0, 1], 10]}, 
        model1[α, β, ν, x0]\ ]\)

model1realization (mod1.m)

model1realization::usage="model1realization[α,β,n,σ] generates a realization of model 1 with gaussian measurement noise with standard deviation σ. This version is based on Matlab code by Daniel Kaplan."

\(model1realization[α_, β_, n_Integer, σ_?Positive]
         := 
With[{
          foo\  = \ RandomArray[NormalDistribution[0, σ], n]}, 
\ 
        model1[α, β, foo, 
          β/\((.00001 + Abs[1 - α])\)]\ ]\)

model2

model2::usage = 
      "model2[α,β,ω,ν,x0] runs model 2 dynamics and \
measurements given the noisy dynamics ω and measurement inputs ν. \
x0 is the initial condition. This version is based on Matlab code by Daniel \
Kaplan."

\(\(model2[α_, β_, ω_?VectorQ, ν_?VectorQ, x0_]
         := 
	ν\  + 
        Rest[FoldList[α\ *\ #1 + β\  + #2\ \ &, x0, ω]]\ 
    \)\)

\(model2[α_, β_, ω_?VectorQ, ν_?VectorQ] := 
	
      With[{x0\  = \ β/\((1 - α)\)}, 
        model2[α, β, ω, ν, x0]\ ]\)

\(model2[α_, β_, ω_?VectorQ] := 
	
      With[{x0\  = \ β/\((1 - α)\), 
          ν = RandomArray[NormalDistribution[0, 1], Length[ω]]}, 
        model2[α, β, ω, ν, x0]\ ]\)

\(model2[α_, β_] := 
	
      With[{x0\  = \ β/\((1 - α)\), 
          ω = \ RandomArray[NormalDistribution[0, 1], 10], 
          ν = RandomArray[NormalDistribution[0, 1], 10]}, 
        model2[α, β, ω, ν, x0]\ ]\)

model2realization (mod2.m)

model2realization::usage="model2realization[α,β,n,sdyn,smeas] generates a realization of model 2, where n = number of points, sdyn = standard deviation of dynamical noise, smeas = standard deviation of measurement noise. This version is based on Matlab code by Daniel Kaplan."

\(model2realization[α_, β_, n_Integer, sdyn_?Positive, 
        smeas_?Positive] := 
	
      With[{ndyn\  = RandomArray[NormalDistribution[0, sdyn], n], 

          nmeas\  = RandomArray[NormalDistribution[0, smeas], n]}, 
\ 
        model2[α, β, ndyn, \ nmeas, \ 
          β/\((.00001 + Abs[1 - α])\)]]\)

crinkleStatistic

crinkleStatistic::usage="crinkleStatistic[y] calculates the crinkle statistic on a vector of de-meaned data:
(x_{t-1}-2*x_t+x_{t+1})^4> / < ( x_t^2 ) ^2
	as proposed by James Theiler.
This version is based on Matlab code by Daniel Kaplan."

\(crinkleStatistic[data : {__?NumericQ}] /; \((Length[data] > 3\ )\) := 

      With[{x\  = data\  - Mean[data]}, 
        With[{x2 = \((Mean[x^2])\)^2}, 

          If[x2\  == \ 0\ , 0, 
				
            With[{d2\  = \ 
                  2*Drop[Rest[x], \(-1\)] - \ Drop[x, \(-2\)]\  - \ 
                    Drop[x, 2]}, 
\ Mean[d2\ ^\ 4]/x2]]\ ]]\)

kaplanScatterPlot

kaplanScatterPlot::usage="kaplanScatterPlot[data,tau,opts]makes a scatter plot of x(t+tau) versus x(t). This version is based on Matlab code by Daniel Kaplan."

\(kaplanScatterPlot[data : {__?NumericQ}, τ_:  1, opts___Rule] /; 
        \((IntegerQ[τ] && Positive[τ] && 
            \((τ < Length[data])\))\) := 

      ListPlot[Transpose[{Drop[data, τ], Drop[data, \(-τ\)]}], 
        PlotRange -> All, AspectRatio -> 1, 
        FilterOptions[ListPlot, opts]]\)

\(kaplanScatterPlot3[data : {__?NumericQ}, τ_:  1, opts___Rule] /; 
        \((IntegerQ[τ] && Positive[τ] && 
            \((2\ τ < Length[data])\))\) := 

      ScatterPlot3D[
        Transpose[{Drop[data, 2\ τ], 
            Drop[Drop[data, τ], \(-τ\)], 
            Drop[data, \(-2\)\ τ]}], PlotRange -> All, 
        FilterOptions[ScatterPlot3D, opts]]\)

nonlinearPredictionError  (prederr.m) (still not working)

Options[nonlinearPredictionError] = {EmbeddingLag -> 1, 
        EmbeddingDimension -> 3}

nonlinearPredictionError::usage="nonlinearPredictionError[x] returns the nonlinear prediction error for time series x. This version is based on Matlab code by Daniel Kaplan."

EmbeddingDimension::usage = 
      "EmbeddingDimension is an option to nonlinearPredictionError with a \
default value of 3."

EmbeddingLag::usage = 
      "EmbeddingLag  is an option to nonlinearPredictionError with a \
default value of 1."

\(nonlinearPredictionError[x_?VectorQ, opts___Rule] := 

      Module[{data2, images, lin, const}, 
		
        With[{lag = 
              \(EmbeddingLag /. {opts}\) /. 
                Options[nonlinearPredictionError], 
				
            emdim = \(EmbeddingDimension /. {opts}\) /. 
                Options[nonlinearPredictionError]}, 
		
          With[{k = If[emdim > 19, emdim + 1, 20]}, 
            
{data2, \ images}\  = \ 
              getImage[\ mHistory[x, emdim, \ lag], \ lag]; 
            
{lin, \ const}\  = \ evalMap[\ data2, \ images, \ k]; 
\ 
            Standardize[images - lin]\ ]]]\)

\(\(nonlinearPredictionError[x_?VectorQ, a_Integer, b_Integer] := 

      nonlinearPredictionError[x, EmbeddingLag -> a, 
        EmbeddingDimension -> b]
\)\)

nearNeighbour (nearneib.m)

The function nearneib.m doesn't appear to do anything, given the content of the original code.  I have implemented the row number (actually column in this case) test and the pattern for the function to match, but there is currently no return value.

nearNeighbour::usage="nearNeighbour[emb, p, k, time, lockout] finds the k nearest neighbors of a point in a set of vectors where emb is the list of vectors (embedded points?), p is the reference point, k is the number of neighbors and lockout is an optional paramater that excludes points that are nearby in time. This version is based on Matlab code by Daniel Kaplan."

nearNeighbour[emb_?MatrixQ,  p_List,  k_,  time_,  lockout_] /; 
        (Dimensions[emb][[2]] == 
            Length[p] )  :=Null

polynomialFitwithCrossValidation  (polycv.m)

polynomialFitwithCrossValidation::usage="polynomialFitwithCrossValidation[pre,post,order] calculates the residuals of a polynomial fit of post vs pre using leave-one-out cross validation.
This version is based on Matlab code by Daniel Kaplan."

polynomialFitwithCrossValidation[pre_?VectorQ,post_?VectorQ,order_Integer]/;
(Length[pre]==Length[post]):=
	With[{npts=Length[pre],inds=Range[Length[pre]],vars=Table[x^n,{n,0,order}]},
(* loop over all points *)
	Table[(Fit[Transpose[{Delete[pre,k],Delete[post,k]}],vars,x]/.x->pre[[k]])-post[[k]],{k,npts}] ]

powerTest

powerTest::usage = 
      "powertest[x,n] calculates Q=<(x_{t+1})^n> normalized by <x^2>^{n/2} \
on a vector 
This version is based on Matlab code by Daniel Kaplan."

\(powerTest[x_?VectorQ, n_] := 
	
      With[{x2\  = \ Mean[x\ *x\ ]^\((n/2)\)}, 
		
        If[x2 == 0, 0, 

          With[{d2 = Rest[x] + Drop[x, \(-1\)], 
              d2a = Rest[x] - Drop[x, \(-1\)]}, {\ Mean[d2\ ^\ n]/x2, 
              Mean[x\ ^\ n]/x2, Mean[d2a\ ^\ n]/x2}\ ]]]\)

squidSinglePulse  (squid1.m)

squidSinglePulse::usage="squidSinglePulse[ amp, dur]
stimulates a model of the Squid Giant Axon with a single pulse of amplitude amp and duration dur,
and returns the voltage sampled at 100 samples per time unit.
This version is based on Matlab code by Daniel Kaplan."

squidSinglePulse[amp_?Positive, dur_?Positive]:=
With[{a = 5 + dur + 30,stims={5,amp,dur} },
(* save /tmp/stims1 stims -ascii *)
(* eval(['! /ssnd/bin/hh ' num2str(a) ' 1 < /tmp/stims1 > /tmp/res.dat']); *)
(* load /tmp/res.dat; *)
TakeColumns[res,{3,3}]]

squidTwoPulse (squid2.m)

squidTwoPulse::usage="squidTwoPulse[ amp1, dur1, interval, amp2, dur2]
does 2-pulse stimulation of the squid model where amp1, dur1 = amplitude, duration  of the first pulse; interval = interval between the start of the two pulses; and amp2, dur2 = amplitude, duration  of the second pulse.
This version is based on Matlab code by Daniel Kaplan."

squidTwoPulse[ amp1_?Positive, dur1_?Positive, interval_?Positive, amp2_?Positive, dur2 _?Positive]:=
	With[{a = 5 + dur1 + dur2 + interval + 30,
stims={{5,amp1,dur1},{5+interval,amp2,dur2}} },
(* save /tmp/stims2 stims -ascii *)
(* eval(['! /ssnd/bin/hh ' num2str(a) ' 1 < /tmp/stims2 > /tmp/res.dat']); *)
(* load /tmp/res.dat; *)
TakeColumns[res,{3,3}]]

timeReversalAssymetryStatistic (timerev.m)

timeReversalAssymetryStatistic::usage = 
      "timeReversalAssymetryStatistic[x,lag
calculates a time reversal asymmetry statistic.
This version is based on Matlab code by Daniel Kaplan."

Here's an faster version that doesn't require loading in LinearAlgebra`MatrixManipulation`

\(timeReversalAssymetryStatistic[x_, lag_:  1] := 
	
      With[{foo = mHistory[x, 3, lag]}, 
		
        With[{a = 
              \(Transpose[foo]\)[[1]], 
            b = \(Transpose[foo]\)[[2]], 
            c = \(Transpose[foo]\)[[3
                ]]}, 
Mean[a*a*b\  - \ b*c*c]]]\)

Here's an alternative version that uses LinearAlgebra`MatrixManipulation`, which looks neater, but is slower in the end.

\(timeReversalAssymetryStatisticSlower[x_, lag_:  1] := 
	
      With[{foo = mHistory[x, 3, lag]}, 
		
        With[{a = Flatten[TakeColumns[foo, {1, 1}]], 
            b = Flatten[TakeColumns[foo, {2, 2}]], 
            c = Flatten[TakeColumns[foo, {3, 3}]]}, 

          Mean[a*a*b\  - \ b*c*c]]]\)

theilerQ

theilerQ::"usage"="theilerQ calculates Q=(x_t +x_{t+1})^3 normalized by x^2^{3/2} where x is vector. This version is based on Matlab code by  Daniel Kaplan."

\(theilerQ[x : {__?NumericQ}] := 

      With[{x2\  = Mean[x\ *\ x]^\((3/2)\)}, 
		
        If[x2 == 0, 0, Mean[\((Rest[x] + Drop[x, \(-1\)])\)^3]/x2]]\)

typicalDistances (disttyp.m)

Clear[typicalDistances]

typicalDistances::usage="typicalDistances[z,percs] calculates typical distances between pre-images where z is some embedded data, and percs are the percentiles to use. It returns distances at the given percentiles. All of this is from a small sample of all pairs of distances. This version is based on Matlab code by  Daniel Kaplan."

typicalDistances[z_,percs_]:=
(* all the elements of percs must be between 0 and 1 *)
Module[{ntargs,dists},With[{Npts = Length[z]},
(* we want up to 10000 points *)
ntargs = Min[Npts, Ceiling[10000/(Npts-1)]];
 dists=Drop[Sort[Take[Flatten[Table[oneDist[z,z[[j]]],{j,ntargs}]],Npts*ntargs]],ntargs];
(* get rid of the zeros corresponding to self-distances *)
With[{𝒹=Length[dists]}, Table[dists[[ Min[ Max[Round[percs[[p]]*(𝒹-1)],0],𝒹-1]+1 ]],{p,1,Length[percs]}] ]]] /;(And @@ (((# >0. )&& (# <1.))& /@ percs))

zeroPredictionError (zeropred.m) (not quite finished?)

zeroPredictionError::usage="zeroPredictionError calculates prediction error from time-series. Calculates mean error of zeroth order prediction of entire series. Mean error is normalized by dividing by error of constant
prediction = mean of series, to get NPE (norm. pred. error). This version is based on Matlab code by  Daniel Kaplan. This Mathematica version has four options with the following defaults: {NumberOfNeighbours->10, MaxHorizon->10, Lag->1, EmbeddingDimension->3}"

Options[zeroPredictionError] = {NumberOfNeighbours -> 10, 
        MaxHorizon -> 10, Lag -> 1, EmbeddingDimension -> 3}

zeroPredictionError[data_?VectorQ, opts___Rule]:=With[\({𝒩 = Length[data], 
			
            k = \(NumberOfNeighbours /. {opts}\) /. 
                Options[zeroPredictionError], 
			
            Hmax = \(MaxHorizon /. {opts}\) /. Options[zeroPredictionError], 
            
			ℒ = 
              \(Lag /. {opts}\) /. Options[zeroPredictionError], 
			
            ℰ = 
              \(EmbeddingDimension /. {opts}\) /. 
                Options[zeroPredictionError], 
			
            huge = 1000000000000000000}\),
Module[{xx, me, End, m1, m2, normE, a, tot,  error,  merror, sorted, 
                index},
			xx = mHistory[data, ℰ, ℒ] ;
  
   me = Mean[
                    Take[Transpose[xx]
                          [[ℰ]], {
                        1 + (ℰ - 1)*ℒ + 
                          2*Hmax, 
                        𝒩 - 
                          (ℰ - 1)*ℒ}]];
   End = 
                  𝒩 - 
                    (ℰ - 1)*ℒ - Hmax;
  
Do[   
       m1 = Max[0, j - 50]; m2 = Min[j + 50, End + 1];
      normE = 
                        Join[Table[
                            (xx[[i, 
                                    Range[ℰ]
                                      ]] - 
                                  xx[[j, 
                                    Range[ℰ]
                                      ]]).
                              (xx[[i, 
                                    Range[ℰ]
                                      ]] - 
                                  xx[[j, 
                                    Range[ℰ]
                                      ]]), {i, m1}], 
                          Table[huge, {m2 - m1}]];
      {sorted, index} = Sort[normE];
		     
                          
  
  
  
       error[[j]] = 
                        Table[xx[[j + H, 
                              ℰ]] - 
                            (Plus @@ 
                                  ( xx
                                      [[Range[H + 1, 
                                        H + k], 
                                      ℰ]])
                                    /k), {H, Hmax}];
       merror[[j]] = 
                        Table[xx[[j + H, 
                              ℰ]] - me, {
                            H, Hmax}];   ,{j, End}];      

  
     npe = 
                  Table[ norm[
                        Normalize[error]
                            [[𝒽]]]/ 
                      Normalize[
                        Transpose[merror]
                            [[𝒽]]], {𝒽, Hmax}];
  
Transpose[{Range[Hmax], 
                    TakeMatrix[ npe, {1, 1}, {Hmax, ℰ}]}];
]]

Maps

Maps are also really easy to do as anonymous functions.

henonMap

henonMap::usage="henonMap[{x1,x2},p] gives a single iteration of the Henon map. Note that the p parameter doesn't appear to do anything and is included only for comparability with other map functions in this notebook.This version is based on Matlab code by Daniel Kaplan."

\(henonMap[{x1_, x2_}, p_] := 
      
{1.4\  + \ 0.3\ *\ x2\  - \ x1*x2, \ x1}\)

quadraticMap

quadraticMap::usage = 
      "quadraticMap[x,p] gives a single iteration of the quadratic map. \
This version is based on Matlab code by Daniel Kaplan. The parameter p is \
optional and has a default value of 4, the chaotic value."

quadraticMap[x_, p_:  4] :=  p * x * (1 - x);

squidMap

squidMap::usage="squidMap[x,p] creates something that looks like the squid map, symmetric around x=0.5, where x is the state
and p is a parameter representing inter-stimulus interval. This version is based on Matlab code by Daniel Kaplan. ALERT! Make sure you use the appropriate values for x and p  to avoid Overflow."

\(squidMap[x_, p_] := 
	
      With[{𝓍 = If[x > 0.5, 1 - x, x]}, 

        If[\ 𝓍\  < \ .693/p, \((\ Exp[𝓍*p] - 1)\)^3, \ 
          1 - .05*𝓍]]\)

tentMap

tentMap::usage="tentMap[x,p] tent map of x at parameter p.
p = 1.99999 is the chaotic one. This version is based on Matlab code by Daniel Kaplan."

tentMap[x_, p_?NumericQ] /; (x > 0  &&  x < 1) := 
      p* If[x  >  0.5, 1 - x, x]

ListPlot[NestList[tentMap[#, 1.8]&, 0.49, 1000], 
        PlotJoined -> True] // Timing

{0.266667 Second,⁃Graphics⁃}

{0.333333 Second,⁃Graphics⁃}

delta-epsilon

getImage

getImage::usage="getImage[data, pred]
finds the scalar images of the points in a time series <pred> time sets in the future
data --- matrix of embedded data (from mHistory)
pred --- look ahead time, default value 1.  It returns data2 --- a new embedded data matrix appropriately trimmed, and images --- the images (at time <pred>) of the points in data2. This is a convenience program to trim an embedding appropriately.
 This version is based on Matlab code by  Daniel Kaplan."

\(getImage[data_?MatrixQ, pred_Integer:  1] := 
      
{\ Drop[data, pred], First\ /@\ Drop[data, \(-pred\)]}\)

\(getImage[data_List, pred_Integer:  1] := 
      
{\ Drop[data, pred], First\ /@\ Drop[data, \(-pred\)]}\)

oneDist

oneDist::usage="oneDist[data,pt] calculates the distance between point pt and each row in matrix data.
This version is based on Matlab code by  Daniel Kaplan."

\(oneDist[data_?MatrixQ, pt_?VectorQ] /; 
        \((\(Dimensions[data]\)[[2]]\ 
             == \ Length[pt])\)\  := 

      Sqrt[Table[
          Plus\ @@\((
              \((data[[i]] - pt)\)^2)
              \), {i, Length[data]}]]\)

\(oneDist[data_?VectorQ, pt_?NumericQ] := 

      Sqrt[Table[
          Plus\ @@\((
              \((data[[i]] - pt)\)^2)
              \), {i, Length[data]}]]\)

deltaepsilonFit

deltaepsilonFit::usage="linear fitting routine for delta-epsilon method.
delta   --distances between pre-images: output by delta-epsilon
epsilon -- distances between images: output by delta-epsilon
maxdelta-- optional - largest delta to consider.
This version is based on Matlab code by  Daniel Kaplan."

deltaepsilonFit[delta_,epsilon_,maxdelta_:(10^100)]:=
(* get rid of points whose pre-images are too far apart *)
Module[{x,y,i,s,sx,sy,sxx,sxy,del},
With[{foo =Transpose[Select[ Transpose[delta],#[[2]] < maxdelta &]]},
		With[{delta2 = delta[[foo]],epsilon2 = epsilon[[foo]]},
(* Sort the pairs into ascending order on delta *)
{x,i}= Sort[delta2];
y = epsilon2[[i]];
s = Length[x];
sx =Plus @@ x;
sy =Plus @@ y;
sxx =Plus @@( x^2);
sxy = Plus @@( x*y);
del = s*sxx - sx*sx;
{(sxx*sy - sx*sxy)/del,(s*sxy - sx*sy)/del} ]]]

deltaEpsilonMethod

deltaEpsilonMethod::usage="deltaEpsilonMethod[z, images, lockout] implements the Delta-epsilon method where z is embedded data as from getImage[], images as from getImage[], and lockout is an optional argument that excludes points closer than the specified value.This version is based on Matlab code by  Daniel Kaplan."

deltaEpsilonMethod[z_,  images_,  lockout_:  0]:=
	Module[{delta, epsilon,  count, kk},
With[{npts  = Length[z]},
			With[{k = Min[npts - 2, 20]},
delta = Flatten[ZeroMatrix[k*npts, 1]];
epsilon = delta;
count = 0;
Do[ (* find the nearest neighbors,eliminating the ones closer in time than 
                              the lockout  *)
  {inds, dists}  = 
                            Transpose[
                              Cases[Transpose[ 
                                  findNeighbour[z,  
                                    z[[i
                                      ]],  k + 10]], 
                                ( Abs[#[[1
                                        ]]  -  i] ) > 
                                     lockout  &] ];
  kk  =  Min[Length[dists], k];
	    
  delta[[count + 1, 
                                  count + kk]]  =  
                                Take[dists, kk];
	    
                                  
\(\ \ 
                              epsilon[[count + 1, 
                                  count + kk]]\  = \ 
                                images
                                    [[Take[inds, kk]
                                      ]]\  - \ 
                                  images[[i
                                    ]]; 
\ \ 
                              \(count\  = \ count + kk; \)\)
					,{i, npts}];
  
{Abs[Take[epsilon, count]], Take[ delta, count]} ]]]

DoDE

DoDE::usage="DoDE[ts,dim,lag,nsurrs] runs deltaeps on some data and surrogates. This version is based on Matlab code by Daniel Kaplan."

\(DoDE[ts_, dim_, lag_, nsurrs_] := 
      
{\ quickDE[ts, dim, lag], 
        Table[quickDE[ampsurr[ts], dim, lag], {nsurrs}]}\)

quickDE

quickDE::usage = 
      "quickDE[ts,dim,lag,nmin] does a quick-and-dirty characterization of \
determinism using de:
ts -- the time series
dim -- the embedding dimension
lag -- the embedding lag
nmin -- optional: number of points to use for delta-eps fitting -- default \
value: 500
This version is based on Matlab code by  Daniel Kaplan."

quickDE[ts_,dim_,lag_,nmin_:500]:=
	Module[{pre,post,perc,deltamax,delta,epsilon},
(* embed the data and get the 'forecasting' pre-image and image *)
{pre,post}= getImage[ mHistory[ts,dim,lag],1];
(* find a 'typical' distance to use for delta-epsilon
 we want to have at least nmin points, so we pick *)
perc = (2*nmin)/(Length[ts])^2;
deltamax = typicalDistances[pre,perc];
{delta,epsilon} = deltaEpsilonMethod[pre,post];
 deltaepsilonFit[delta,epsilon,deltamax] ]

plotDE

plotDE::usage="plotDE[delta,epsilon,opts] makes a cumulative plot of epsilon vs delta that is useful for picking length scales and for comparing to surrogate data. This version is based on Matlab code by  Daniel Kaplan."

plotDE[delta_, epsilon_, opts___Rule]:=
  
\(With[{{d, inds}\  = Transpose[Sort[Transpose[delta]]]}, 
		
          With[{ee\  = \ 
                NestList[Plus, 
                    First[epsilon[[inds
                        ]]], 
                    Rest[epsilon[[inds
                        ]]]]/Range[Length[inds]]}, 

            ListPlot[d, ee, PlotJoined -> True, 
              FilterOptions[ListPlot, opts]]\ ]]\);