1
## Copyright (C) 1996-2015 Piotr Held
3
## This file is part of Octave.
5
## Octave is free software; you can redistribute it and/or
6
## modify it under the terms of the GNU General Public
7
## License as published by the Free Software Foundation;
8
## either version 3 of the License, or (at your option) any
11
## Octave is distributed in the hope that it will be useful,
12
## but WITHOUT ANY WARRANTY; without even the implied
13
## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14
## PURPOSE. See the GNU General Public License for more
17
## You should have received a copy of the GNU General Public
18
## License along with Octave; see the file COPYING. If not,
19
## see <http://www.gnu.org/licenses/>.
22
## @deftypefn{Function File} {[@var{surro_data}, @var{pars}] =} surrogates (@var{S})
23
## @deftypefnx{Function File} {[@var{surro_data}, @var{pars}] =} surrogates (@var{S}, @var{paramName}, @var{paramValue}, @dots{})
25
## Generates multivariate surrogate data (implements the iterative Fourier scheme).
26
## Surrogate data is generated from a dataset with the aim of testing whether the
27
## dataset was generated by a given process (null hypothesis). The Fourier scheme
28
## assumes that the dataset is the output of a Gaussian linear stochastic process.
29
## Surrogate data is generally used to test the null hypothesis.
35
## This function always assumes that each time series is along the longer
36
## dimension of matrix @var{S}. It also assumes that every dimension
37
## (counting along the shorter dimension) of @var{S} is considered a
38
## component of the time series. It's length must be factorizable by only
39
## 2, 3 and 5. If not the largest submatrix that fulfills this requirement will be
40
## used. The function @code{endtoend} can be used to determine what is the best
41
## submatrix for the data and then sending only that submatrix to this program.
42
## Padding with zeros is @strong{not} and option.
45
## @strong {Parameters}
49
## Sets the number of surrogates to be calculated. Determines the form of
50
## the output (see Output section) [default = 1].
52
## The maximum number of permutations. Value '0' yields random permutations or if
53
## switch @var{exact} is set an unrescaled FFT surrogate. Value '1' is a surrogate
54
## close to the result of the AAFT procedure, but not quite the same. Value '-1'
55
## means the program will perform iterations until there is no change between them
58
## Set the seed for the random generator [default = use default seed].
65
## This switch makes the spectrum of the output exact rather than a distribution.
72
## If parameter @code{n == 1} then this is a matrix that holds the surrogate data.
73
## If parameter @code{n > 1} then it is @var{n} x 1 cell array of matrixes with
74
## the data. In both cases the matrixes themselves are alligned with the input.
76
## This is a matrix of size @var{n} x 2 (if the input components were column
77
## vectors, otherwise transposed). The first column contains the number of
78
## iteration it took to generate the @var{i}-th surrogate, whereas the second
79
## column is the relative discrepency for the @var{i}-th surrogate.
82
## @seealso{demo surrogates, endtoend}
84
## @strong{Algorithms}
86
## The algorithms for this functions have been taken from the TISEAN package.
89
## Author: Piotr Held <pjheld@gmail.com>.
90
## This function is based on surrogates of TISEAN 3.0.1
91
## https://github.com/heggus/Tisean"
93
function [surro_data, pars] = surrogates (S,varargin)
95
if (nargin < 1 || nargout > 2)
99
if ((ismatrix (S) == false) || (isreal(S) == false) || ...
100
(isreal(S) == false))
101
error ('Octave:invalid-input-arg', "S must be a realmatrix");
106
max_iterations = -1; # means until no change
111
p.FunctionName = "surrogates";
113
isPositiveIntScalar = @(x) isreal(x) && isscalar (x) ...
114
&& (x > 0) && (x-round(x) == 0);
115
isValidIterationValue = @(x) isPositiveIntScalar (x) ...
116
|| (isscalar(x) && ((x == 0) || (x == -1)));
117
isNonNegativeScalar = @(x) isreal(x) && isscalar (x) ...
118
&& ((x > 0) || (x == 0));
120
p.addParamValue ("n", nsur, isPositiveIntScalar);
121
p.addParamValue ("i", max_iterations, isValidIterationValue);
122
p.addParamValue ("seed", seed, isNonNegativeScalar);
123
p.addSwitch ("exact");
125
p.parse (varargin{:});
129
max_iterations = p.Results.i;
130
seed = p.Results.seed;
131
ispec = p.Results.exact;
133
# Check if nsur is not large
135
warning ("Octave:tisean", ["Parameter 'n' is larger than 1000, this ", ...
136
"function migth execute a long time and take ",...
137
"up a lot of memory"])
140
# Correct S to always have more rows than columns
142
if (rows (S) < columns (S))
147
# Check if the length of 'S' can be factorized by only 2, 3 and 5
148
original_length = length (S);
149
while (max (factor (length (S))) > 5)
153
if (original_length > length (S))
154
warning ("Octave:tisean",...
155
["The length of 'S' was not factorizable by 2, 3 and 5. ", ...
156
"Using only first %d so that it's length would fulfill this ",...
157
"requirement"], length (S));
161
[surro_data, pars] = __surrogates__ (S, nsur, max_iterations, ispec, seed);
163
# If the input was transposed allign output with it
166
surro_data = cellfun (@transpose, surro_data, 'UniformOutput', false);
169
# If surro_data is a single cell, then change it to be a matrix
170
if (isequal (size (surro_data), [1,1]))
171
surro_data = surro_data {1};
177
%% 'x' will be a stationary Gaussian linear stochastic process
178
%! x = zeros (2000,1);
180
%! x(i) = 0.7*x(i-1) + (-6 + sum (rand ([size(1), 12]), 3));
183
%! # 'spike' is the process above measured s_n (x_n) = x_n^3.
192
%! plot (surrogates(spike),'b');
194
%! title ("surrogates")
195
%!###############################################################
198
%! [s,p] = surrogates (henon (1000).', 'i', 50, 'n', 15);
201
%% Check if parameter i (max iterations works properly)
203
%! expected(1:15) = 50;
204
%! assert (p(:,1), expected.')
206
%% Check if the relative discrepancy remains similar
207
%!assert (std (p(:,2)), 0, 5e-4)
209
%% Check if cell was properly created and transposed
210
%!assert(iscell (s) && isequal (size (s), [15, 1]))
211
%!assert(rows (s{1}) < columns (s{1}))
213
%% Check if shortening data is discovered
214
%% Warnings are promoted to errors to avoid further computation
215
%!error <factorizable> warning ("error", "Octave:tisean"); surrogates (henon (11));
217
%% Check if shortening the data works as expected
220
%! expected = surrogates (henon (100), 'seed', 0.25);
221
%! result = surrogates (henon (102), 'seed', 0.25);
222
%! assert (result, expected)
224
%% Check if checking parameter 'n' works as expected
225
%!error <long> warning ("error", "Octave:tisean"); surrogates ([1 2], 'n', 1001);
added
removed
inst/surrogates.m
