-
Notifications
You must be signed in to change notification settings - Fork 2
/
corrint.m
executable file
·309 lines (296 loc) · 8.89 KB
/
corrint.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
function varargout=corrint(varargin)
%
% [y1,y2,y3]=corrint(x,embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode,findScaling)
%
% Correlation integral analysis of a time series. Based on:
%
% [1] Kaplan, Daniel, and Leon Glass. Understanding nonlinear dynamics. Vol. 19. Springer, 1995.
% [2] Kantz, Holger, and Thomas Schreiber. Nonlinear time series analysis. Cambridge university press, 2004.
%
% Required input parameter:
% x
% Nx1 matrix (doubles) of time series to be analyzed.
%
% Optional Parameters are:
%
% embeddedDim
% 1x1 Integer specifying the embedded dimension size to use (default
% =2).
%
% timeLag
% 1x1 Integer specifying the minimum time lag distance (in samples) of the point to
% be estimated. Default is 2. If timeLag=-1 the timeLag is estimated
% from the first zero-crossing point of the autocorrelation of x.
%
% timeStep
% 1x1 Integer specifying time lag distance (in samples) within
% each point used in the embeddedDimm vector. For example, if embeddedDim
% is 3 and timeStep =2, then the embedded dimension vector will consists of
% 3 samples separated by 2 samples each, covering a window of size of 7 samples.
%
% distanceThreshold
% 1x1 double specifying the distance threshold between embedded
% points. The points who's distance is less than distanceThreshold are considered in
% the same neighborhood and used for either prediction, recurrence, or the
% estimation of the embedded dimension.
%
% neighboorSize
% 1x1 Integer specifying the number of neighbors to be used for
% prediction and smoothing (see 'estimationMode' parameter).
%
% estimationMode
% String specifying what analysis type to be done in the time series.
% Options are:
% 'recurrence' -Calculates recurrence data to be used
% in for recurrence plots (default).
% 'dimension' -Generates statistics for the estimation of the correlation dimension
% of the time series and it's scaling
% regions.
% 'prediction' -Predicts second half of the time
% series using the first half as a model
% and neighboorSize nearest points.
% 'smooth' -Predicts all point of the times
% series using all other points as a
% model and neighboorSize nearest
% points.
%
% findScaling
% 1x1 Boolean flag to be passed when using 'dimension' mode. If set to
% true, the scaling region will be searched automatically, using
% r1=std(x)/4 and r2 -> C(r1)/C(r2) ~ 5. Default value is false.
%
% The output returned by CORRINT is dependent on the 'estimationMode'
% parameter, so that the description of the output below is broken down into the
% different possible options for the 'estimationMode' parameter.
%
% Output Parameters - 'recurrence' mode
% y1
% Lx1 Vector of integers for state i.
% y2
% Lx1 Vector of integers for state state j that is a neighbor of state i (first column).
%
% Output Parameters - 'dimension' mode
% y1
% Lx1 Vector of doubles of log(distanceThreshold).
% y2
% Lx1 Vector of doubles for log(neighborhood size) given the distanceThreshold used in column 1.
% y3
% 1x1 double. Optional, estimated slope of y1 and y2
%
% Output Parameters - 'prediction' mode
% y1
% Lx1 Vector of doubles of estimated second half of the time series.
% y2
% Lx1 Vector of doubles for original second half of the time series.
% y3
% 1x1 double. Optional, variance of the prediction error divided by variance of the second half of the time series.
%
%
% Output Parameters - 'smooth' mode
% y1
% Lx1 Vector of doubles of smoothed the time series.
% y2
% Lx1 Vector of doubles for original time series.
% y3
% 1x1 double. Optional, variance of the prediction error divided by variance of the time series.
%
%
% %%% Beging Example %%%
%
% N=500; %Number of points for each process
% model_names={'linearModel','nonlinearModel'};
%
% %Linear Auto Regressive model with measurement noise
% linearModel=zeros(N,1);
% x=77;
% linearModel(1)=x;
% for n=2:N
% x=4 + 0.95*x;
% linearModel(n)= x + randn(1)*2;
% end
%
% %Non-linear model of dimension ~ 3.9
% nonlinearModel=zeros(N,1);
% x=0.2;y=0.2;z=0.2;v=0.2;model_five(1)=x;
% for n=2:N
% m=0.4 - 6/(1+ x^2 + y^2);
% xold=x;yold=y;zold=z;vold=v;
% x= 1 + 0.7*(xold*cos(m)-yold*sin(m)) + 0.2*zold;
% y=0.7*(xold*sin(m) + yold*cos(m));
% z=1.4 + 0.3*vold - zold^2;
% v=zold;
% nonlinearModel(n)= x + 0.3*z + randn(1)*0.05;
% end
%
% %Plot time series
% figure(1)
% for i=1:2
% subplot(2,1,i)
% eval(['plot(' model_names{i} ');legend(''' model_names{i} ''')'])
% title('Time Plot');xlabel('time')
% end
%
% %Plot cross correlation
% figure(2)
% for i=1:2
% subplot(2,1,i)
% eval(['x=' model_names{i} ';'])
% R=xcorr(x-mean(x),'coeff');
% plot(R(round(N):end))
% eval(['legend(''' model_names{i} ''')'])
% title('Autocorelation'); xlabel('lag')
% end
%
% %Plot Phase Plots
% figure(3)
% for i=1:2
% subplot(2,1,i)
% eval(['x=' model_names{i} ';'])
% scatter(x(1:end-1),x(2:end))
% eval(['legend(''' model_names{i} ''')'])
% title('Phase Plot');xlabel('x(t)');ylabel('x(t+1)')
% end
%
% %Plot prediction errors vs surrogate
% timeLag=1;
% timeStep=1;
% distanceThreshold=[];
% embeddedDim=4;
% estimationMode='smooth';
% figure(4)
% K=[1:20 25 30 50 70 100];
% D=length(K);
% surrN=10;
% for i=1:2
% eval(['x=' model_names{i} ';'])
% err=zeros(D,1)+NaN;
% surr_data=zeros(D,surrN);
% SURR=surrogate(x,surrN);
% for d=1:D;
% neighboorSize=K(d);
% [y1,y2,y3]=corrint(x,embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode);
% err(d)=y3;
% for s=1:surrN
% [y1,y2,y3]=corrint(SURR(:,s),embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode);
% surr_data(d,s)=y3;
% end
% end
% subplot(2,1,i)
% plot(K,err,'o-');hold on
% errorbar(K,mean(surr_data,2),var(surr_data,[],2)./sqrt(10),'r')
% eval(['legend(''' model_names{i} ''',''surrogate'')'])
% xlabel('Embedded Dimension')
% ylabel('err/var')
%
% end
%
% %%% End Example %%%
%
% Written by Ikaro Silva, 20134
% Last Modified: November 23, 2014
% Version 1.0
%
% Since 0.9.8
%
%
% See also SURROGATE, DFA, MSENTROPY
%endOfHelp
persistent javaWfdbExec config
if(isempty(javaWfdbExec))
[javaWfdbExec,config]=getWfdbClass('corrint');
end
%Set default pararamter values
inputs={'x','embeddedDim','timeLag','timeStep','distanceThreshold','neighboorSize','estimationMode','findScaling'};
outputs={'y1','y2','y3'};
embeddedDim=[];
timeLag=[];
timeStep=[];
distanceThreshold=[];
neighboorSize=[];
estimationMode='recurrence';
findScaling=0;
wfdb_argument={};
y1=[];
y2=[];
y3=[];
for n=1:nargin
if(~isempty(varargin{n}))
eval([inputs{n} '=varargin{n};'])
end
end
if(~isempty(embeddedDim))
wfdb_argument{end+1}='-d';
wfdb_argument{end+1}=num2str(embeddedDim);
end
if(~isempty(timeLag))
wfdb_argument{end+1}='-t';
wfdb_argument{end+1}=num2str(timeLag);
end
if(~isempty(timeStep))
wfdb_argument{end+1}='-s';
wfdb_argument{end+1}=num2str(timeStep);
end
if(~isempty(distanceThreshold))
wfdb_argument{end+1}='-r';
wfdb_argument{end+1}=num2str(distanceThreshold);
end
if(~isempty(neighboorSize))
wfdb_argument{end+1}='-n';
wfdb_argument{end+1}=num2str(neighboorSize);
end
switch estimationMode
case 'recurrence'
wfdb_argument{end+1}='-p';
case 'dimension'
wfdb_argument{end+1}='-v';
y3=' ';
if(findScaling)
wfdb_argument{end+1}='-a';
end
case 'prediction'
wfdb_argument{end+1}='-P';
case 'smooth'
wfdb_argument{end+1}='-S';
y3=' ';
otherwise
error(['Unkown estimation mode: ' estimationMode])
end
javaWfdbExec.setArguments(wfdb_argument);
if(config.inOctave)
x=cellstr(num2str(x));
x=java2mat(javaWfdbExec.execWithStandardInput(x));
Nx=x.size;
out=cell(Nx,1);
for n=1:Nx
out{n}=x.get(n-1);
end
else
out=cell(javaWfdbExec.execWithStandardInput(x).toArray);
end
M=length(out);
if(~isempty(y3))
y3=out{end};
out(end)=[];
M=M-1;
if(strcmp(estimationMode,'smooth'))
tmp=y3;
sep=regexp(tmp,'\s');
y3=str2num(tmp(sep(end):end));
end
end
if(~isempty(strfind(out{1},'Possibly')))
warning(out{1})
out(1)=[];
M=M-1;
end
y1=zeros(M,1)+NaN;
y2=zeros(M,1)+NaN;
for m=1:M
str=out{m};
sep=regexp(str,'\s');
y1(m)=str2num(str(1:sep));
y2(m)=str2num(str(sep(1):end));
end
for n=1:nargout
eval(['varargout{n}=y' num2str(n) ';'])
end