From 5c6feb30b0f3d251e9f712d0bffaefe3cd476ca3 Mon Sep 17 00:00:00 2001
From: "Nes, Egbert van" <egbert.vannes@wur.nl>
Date: Fri, 26 Mar 2021 09:13:51 +0000
Subject: [PATCH] Update adjustplot.m, ebisuzaki.m, fishercombine.m,
generic_ews.m, generic_ews_sens.m, run_regime.m, violinplot.m, README.md
files
---
README.md | 2 +-
adjustplot.m | 383 +++++++++++++++++
ebisuzaki.m | 68 +++
fishercombine.m | 20 +
generic_ews.m | 1008 ++++++++++++++++++++++++++++++++++++++++++++
generic_ews_sens.m | 454 ++++++++++++++++++++
run_regime.m | 81 ++++
violinplot.m | 110 +++++
8 files changed, 2125 insertions(+), 1 deletion(-)
create mode 100644 adjustplot.m
create mode 100644 ebisuzaki.m
create mode 100644 fishercombine.m
create mode 100644 generic_ews.m
create mode 100644 generic_ews_sens.m
create mode 100644 run_regime.m
create mode 100644 violinplot.m
diff --git a/README.md b/README.md
index f50402b..b67d4e7 100644
--- a/README.md
+++ b/README.md
@@ -19,7 +19,7 @@ ebisuzaki.m - generate null models according to the Ebisuzaki (1997)
fishercombine.m - combine independent p values according to the Fisher method
generic_ews.m - the rolling window method for early warning signals (Dakos et al. 2008)
generic_ews_sens.m - sensitivity analysis of generic_ews
-run_regime.m - the main figures are created with this code, usage:
+run_regime.m - the main figures of Scheffer et al. 2021 are created with this code, usage:
run_regime(data,regime), where data should be a MATLAB table with
(at least) two fields:
trees
diff --git a/adjustplot.m b/adjustplot.m
new file mode 100644
index 0000000..1546515
--- /dev/null
+++ b/adjustplot.m
@@ -0,0 +1,383 @@
+function adjustplot(hfig, comm, varargin)
+ %handy to adjust subplot positions or default appearance
+ %adjustplot(gcf,'-d'): set defaults
+ %adjustplot(gcf,'-1',orientation): make axes the same of subplots in a orientation
+ %('x','y','b' (both)
+ %adjustplot(gcf,'-t','tr') add text right (tr), text left (tl), bottom
+ %left (bl), bottom right (bl)
+ %adjustplot(gcf,'-r',orientation): remove axes of a subplot in a
+ %orientation 't'=remove titles only, 'tx'=remove x titles, 'ty' remove
+ %y titles
+ %adjustplot(gcf,'-p',spacing): shift subplots closer together (spacing = [spacingx
+ %spacingy])
+ %adjustplot(gcf,'-b') remove ticks from the axes (does not
+ %automatically update with scaling of the figure)
+ if nargin < 1
+ hfig = gcf;
+ elseif ~ishandle(hfig)
+ if nargin == 2
+ varargin = [{comm} varargin];
+ comm = hfig;
+ hfig = gcf;
+ elseif nargin == 1
+ varargin{1} = [];
+ comm = hfig;
+ hfig = gcf;
+ end
+ end
+ % if nargin < 3
+ % varargin{1} = [];
+ % end
+ %if nargin < 2
+ % comm = '-d';
+ %end
+ switch comm
+ case '-b' %box without ticks, does not work if Inf is in xlim or ylim
+ %not updated (could be done with listener)
+ %ax=get(hfig,'currentaxes');
+ axs = get(hfig, 'children');
+ drawnow()
+ for i = 1:length(axs)
+ if strcmp(get(axs(i), 'Type'), 'axes')
+ yylim = get(axs(i), 'ylim');
+ xxlim = get(axs(i), 'xlim');
+
+ if strcmp(get(axs(i), 'YDir'), 'reverse')
+ iy = [2 1];
+ else
+ iy = [1 2];
+ end
+ if strcmp(get(axs(i), 'XDir'), 'reverse')
+ ix = [2 1];
+ else
+ ix = [1 2];
+ end
+ h = findobj(axs(i), 'tag', 'plot-box');
+ if isempty(h)
+ %first use only!! (else listener causes problem)
+ set(axs(i), 'box', 'off')
+ set(axs(i), 'ylim', yylim);
+ set(axs(i), 'xlim', xxlim);
+ oldhold = ~strcmp(get(axs(i), 'nextplot'), 'add');
+ hold(axs(i), 'on');
+ h = plot(axs(i), [xxlim(ix(1)), xxlim(ix(2)), xxlim(ix(2))], [yylim(iy(2)), yylim(iy(2)), yylim(iy(1))], 'k-', 'tag', 'plot-box', 'LineWidth', get(axs(i), 'LineWidth'));
+ set(get(get(h, 'Annotation'), 'LegendInformation'), 'IconDisplayStyle', 'off');
+ if ~oldhold
+ hold(axs(i), 'off');
+ end
+ addlistener(axs(i), 'YLim', 'PostSet', @(src,evnt)adjustplot(hfig,'-b'));
+ addlistener(axs(i), 'XLim', 'PostSet', @(src,evnt)adjustplot(hfig,'-b'));
+ else
+ set(h, 'xdata', [xxlim(ix(1)), xxlim(ix(2)), xxlim(ix(2))], 'ydata', [yylim(iy(2)), yylim(iy(2)), yylim(iy(1))]);
+ end
+ end
+ end
+ return;
+ case '-d'
+ if nargin < 3
+ apen = struct('pen', '-', 'pen2', '-', 'markersize', 5, 'markerfacecolor', 'none', 'cycle', 0, ...
+ 'i', 1, 'color', [0 0 1], 'color2', [0 0 1], 'drawcolor', [0.7500 0 0], 'colormap', [0 1 1; 1 0 1], ...
+ 'linewidth', 1, 'fontsize', 16, 'tickdir', 'out', 'box', 'on');
+ else
+ apen = varargin{1};
+ end
+
+ try
+ hax = get(hfig, 'CurrentAxes');
+ catch
+ hax = findobj(hfig, 'type', 'axes');
+ end
+ if isempty(hax)
+ hax = gca;
+ end
+
+ Cbar = findobj(hfig, 'tag', 'Colorbar');
+ if ~isempty(Cbar)
+ set(Cbar, 'FontSize', apen.fontsize);
+ end
+ set(hax, 'FontSize', apen.fontsize);
+ set(hax, 'LineWidth', apen.linewidth);
+ set(hax, 'TickDir', 'out');
+ set(hax, 'TickLength', [0.015 0.025])
+ if isempty(findobj(hfig, 'tag', 'plot-box'))
+ set(hax, 'Box', apen.box);
+ else
+ adjustplot(hfig, '-b');
+ end
+ %set(H, 'TickDir', apen.tickdir);
+ % set(hax, 'Box', apen.box);
+ t = findobj(hax, 'type', 'text');
+ if ~isempty(t)
+ if ~isfield(apen, 'fontsize')
+ apen.fontsize = 14;
+ end
+
+ set(t, 'fontsize', apen.fontsize);
+ end
+
+ t1 = [findobj(hfig, 'tag', 'conteq-text'); findobj(hfig, 'tag', 'contbif-text'); findobj(hfig, 'tag', 'label-text')];
+ set(t1, 'fontsize', 10);
+ if isempty(Cbar) || ~any(hax == Cbar)
+ View = get(hax, 'view');
+ if (apen.fontsize >= 14) && min(View == [0 90])
+ set(hax, 'units', 'normalized');
+ if ~isempty(Cbar)
+ P = [0.1500 0.1900 0.6455 0.7350];
+ else
+ P = [0.1300 0.1100 0.7750 0.8150];
+ end
+
+ pos1 = get(hax, 'position');
+ if (pos1(3) > 0.6) && (pos1(4) > 0.6)
+ set(hax, 'position', transform(P, apen.fontsize));
+ if ~isempty(Cbar)
+ P = [0.831375 0.13 0.058125 0.795];
+ set(Cbar, 'position', transform(P, apen.fontsize))
+ end
+
+ end
+
+ end
+ end
+ case '-t'
+ if length(varargin) == 2
+ atext = varargin{2};
+ pos = varargin{1};
+ else
+ atext = varargin{1};
+ pos = 'tl';
+ end
+ hax = get(hfig, 'CurrentAxes');
+ if ischar(pos)
+ switch pos
+ case 'tl'
+ pos = [0.05 0.95];
+ case 'tr'
+ pos = [0.95 0.95];
+ case 'bl'
+ pos = [0.05 0.05];
+ case 'br'
+ pos = [0.95 0.05];
+ end
+ end
+ xlims = get(hax, 'xlim');
+ ylims = get(hax, 'ylim');
+ text(hax, xlims(1) + pos(1) * (xlims(2) - xlims(1)), ylims(1) + pos(2) * (ylims(2) - ylims(1)), atext);
+ case '-l'
+ hs = getsubplotgrid(hfig);
+ if isempty(varargin{1})
+ orientation = 'b';
+ else
+ orientation = varargin{1};
+ end
+ [rows, cols] = size(hs);
+ if strncmpi(orientation, 'b', 1) || strncmpi(orientation, 'x', 1)
+ for i = 1:rows - 1
+ for j = 1:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ adjustticklabels(hs(i, j), 'X')
+ end
+ end
+ end
+ end
+ if strncmpi(orientation, 'b', 1) || strncmpi(orientation, 'y', 1)
+ for i = 1:rows
+ for j = 1:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ adjustticklabels(hs(i, j), 'Y')
+ end
+ end
+ end
+ end
+ case '-r'
+ if isempty(varargin{1})
+ orientation = 'b';
+ else
+ orientation = varargin{1};
+ end
+ hs = getsubplotgrid(hfig);
+ if numel(hs) < 2
+ disp('Too few axes');
+ return;
+ end
+ removeticklabel = true;
+ if strncmpi(orientation, 't', 1)
+ removeticklabel = false;
+ if strfind(orientation, 'x')
+ orientation = 'x';
+ elseif strfind(orientation, 'y')
+ orientation = 'y';
+ else
+ orientation = 'b';
+ end
+ end
+ [rows, cols] = size(hs);
+ if strncmpi(orientation, 'b', 1) || strncmpi(orientation, 'x', 1)
+ for i = 1:rows - 1
+ for j = 1:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ set(get(hs(i, j), 'xlabel'), 'string', '')
+ if removeticklabel
+ set(hs(i, j), 'XTickLabel', []);
+ end
+ end
+ end
+ end
+ end
+ if strncmpi(orientation, 'b', 1) || strncmpi(orientation, 'y', 1)
+ for i = 1:rows
+ for j = 2:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ set(get(hs(i, j), 'ylabel'), 'string', '')
+ if removeticklabel
+ set(hs(i, j), 'YTickLabel', []);
+ end
+ end
+ end
+ end
+ end
+ ndx = hs ~= 0;
+ set(hs(ndx), 'fontsize', 12)
+ aa = get(hs(ndx), 'xlabel');
+ if iscell(aa)
+ aa = [aa{:}];
+ end
+ set(aa, 'fontsize', 12);
+ aa = get(hs(ndx), 'ylabel');
+ if iscell(aa)
+ aa = [aa{:}];
+ end
+ set(aa, 'fontsize', 12);
+ case '-p'
+ if isempty(varargin{1})
+ spacing = 0.02;
+ else
+ spacing = varargin{1};
+ if ischar(spacing)
+ spacing = str2double(spacing);
+ end
+ end
+ set(hfig, 'PaperPositionMode', 'auto')
+ set(hfig, 'units', 'normalized');
+ figpos = get(hfig, 'position');
+ vert = 0.5469; %normal size of figure
+ hor = 0.41;
+ [hs, colwidths] = getsubplotgrid(hfig);
+ if isempty(hs)
+ disp('No axes to adjust');
+ return;
+ elseif length(hs) == 1
+ return;
+ end
+
+ if length(spacing) == 1
+ spacing = spacing + zeros(2, 1);
+ end
+ [rows, cols] = size(hs);
+ if cols <= 4 && rows <= 4
+ newfigpos = [figpos(1:2) hor * cols / 2 vert * rows / 2];
+ else
+ rat = rows / cols;
+ if rat < 1
+ newfigpos = [figpos(1:2) hor vert * rat];
+ else
+ newfigpos = [figpos(1:2) hor / rat vert];
+ end
+ end
+ if newfigpos(4) > 1
+ newfigpos(2) = - 0.2;
+ end
+ set(hfig, 'position', newfigpos);
+ for i = 1:rows
+ for j = 1:cols
+ h = findobj(hfig, 'tag', sprintf('subplot(%d,%d)', i, j));
+ if any(ishandle(h)) && any(h ~= 0)
+ set(h, 'position', subplotpos(rows, cols, j, i, spacing, colwidths(i, j)));
+ end
+ end
+ end
+ movegui(gcf, 'center')
+ otherwise
+ error('adjustsubplots:unknown', 'Unknown command');
+ end
+end
+
+function [hs, colwidths] = getsubplotgrid(h)
+ ch = get(h, 'children');
+ hs = [];
+ colwidths = [];
+ if ~isempty(ch)
+ tags = get(ch, 'tag');
+ ch = ch(~strcmp(tags, 'legend'));
+ types = get(ch, 'type');
+ if ~iscell(types)
+ types = {types};
+ end
+ poss = get(ch, 'position');
+ if ~iscell(poss)
+ poss = {poss};
+ end
+ ipos = zeros(length(ch), 4);
+ for i = length(ch): - 1:1
+ if strcmp(types{i}, 'axes')
+ ipos(i, :) = poss{i};
+ else
+ ipos(i, :) = [];
+ ch(i) = [];
+ types(i) = [];
+ end
+ end
+ %we actually need to account for the spacing between the columns, a
+ %bit complex?
+ colwidth = floor(ipos(:, 3) ./ min(ipos(:, 3)) + 0.01); %the columns each figure spans
+ colpos = sort(unique(sort(ipos(:, 1))), 'ascend');
+ rowpos = sort(unique(sort(ipos(:, 2))), 'descend');
+ hs = zeros(length(rowpos), length(colpos));
+ colwidth(colwidth > length(colpos)) = length(colpos);
+ colwidths = zeros(length(rowpos), length(colpos));
+ for i = 1:length(ch)
+ if strcmp(types{i}, 'axes')
+ arow = find(rowpos == ipos(i, 2));
+ acol = find(colpos == ipos(i, 1));
+ hs(arow, acol) = ch(i);
+ colwidths(arow, acol) = colwidth(i);
+ set(ch(i), 'tag', sprintf('subplot(%d,%d)', arow, acol));
+ end
+ end
+ end
+end
+
+function adjustticklabels(hax, orient)
+ axis(hax, 'tight')
+ newtick = get(hax, [orient 'tick']);
+ tickdiff = (newtick(2) - newtick(1));
+ newlim = [newtick(1) - tickdiff newtick(end) + tickdiff];
+ axis(hax, 'manual')
+ set(hax, [orient 'lim'], newlim);
+ set(hax, [orient 'tick'], [newtick(1) - tickdiff newtick]);
+end
+function P = transform(P, fontsize)
+ %P = [left, bottom, width, height]
+ %where left and bottom define the distance from the lower-left corner of the screen
+ if fontsize <= 18
+ P(1) = P(1) + 0.01;
+ P(2) = P(2) + 0.03;
+ P(3) = P(3) - 0.02;
+ P(4) = P(4) - 0.04;
+ elseif fontsize <= 22
+ P(2) = P(2) + 0.04;
+ P(4) = P(4) - 0.04;
+ elseif fontsize <= 28
+ P(1) = P(1) + 0.02;
+ P(3) = P(3) - 0.02;
+ P(2) = P(2) + 0.08;
+ P(4) = P(4) - 0.08;
+ else
+ P(1) = P(1) + 0.08;
+ P(3) = P(3) - 0.08;
+ P(2) = P(2) + 0.16;
+ P(4) = P(4) - 0.16;
+ end
+
+ return
+end
diff --git a/ebisuzaki.m b/ebisuzaki.m
new file mode 100644
index 0000000..97ffb93
--- /dev/null
+++ b/ebisuzaki.m
@@ -0,0 +1,68 @@
+% X=ebisuzaki(x,nsim,value);
+%
+% This function creates 'nsim' random time series that have the same power
+% spectrum as the original time series 'x' but random phases. The power
+% spectrum of 'x' and 'X' is the same (see abs(fft(x)) and abs(fft(X));
+% The 'X' time series has also the same variance as 'x'. (see Ebisuzaki W. (1997) ;
+% A method to estimate the statistical significance of a correlation when
+% the data are serially correlated. J Climate, 10, 2147-2153).
+%
+% Input
+% 'x' : vector of real number containing the time series to match
+% 'nsim' : integer number giving the number of simulation
+% 'value' : real number to initiate the random sequence (if > 0, the seeds
+% is initiated to the value; otherwise, it is changed from the clock)
+%
+% Output
+% 'X' : matrix of real number (the same number of rows as length of 'x' and
+% 'nsim' columns.
+%
+% fixed by Hao Ye, 2012
+%
+% original version by
+% Vincent MORON
+% March 2002
+function X = ebisuzaki(x,nsim,value)
+if nargin<2
+ nsim=1;
+end
+if nargin<3
+ value=-1;
+end
+
+if(value>0)
+ rand('seed',value);
+else
+ rand('seed',sum(100*clock));
+end
+if any(isnan(x))
+ error('stats:ebisizaki','This procedure cannot handle NaNs');
+end
+
+n = length(x);
+n2 = floor(n/2);
+x = x(:);
+a = fft(x); % discrete fourier transform of the original
+amplitudes = abs(a); % power of the original spectrum
+amplitudes(1) = 0;
+
+for i = 1:nsim
+ if (mod(n, 2) == 0) % series has even length
+ thetas = 2*pi * rand(n2-1,1);
+ angles = [0; thetas(:); 0; flipud(-thetas(:))];
+ recf = amplitudes .* exp(1i*angles); % adjust the phases
+ recf(n2,1) = sqrt(2) * recf(n2,1) * cos(rand()*2*pi); % adjust n/2 power
+ else % series has odd length
+ thetas = 2*pi * rand(n2, 1);
+ angles = [0; thetas(:); flipud(-thetas(:))];
+ recf = amplitudes .* exp(1i*angles); % adjust the phases
+ end
+ X(:,i) = real(ifft(recf));
+end
+
+% adjust variance of the surrogate time series to match the original
+X = X * diag(std(x)./std(X));
+
+% adjust mean of the surrogate time series to match the original
+% X = X - repmat(mean(X), n, 1);
+% X = X + mean(x);
diff --git a/fishercombine.m b/fishercombine.m
new file mode 100644
index 0000000..ff6dd06
--- /dev/null
+++ b/fishercombine.m
@@ -0,0 +1,20 @@
+function [res]=fishercombine(ps)
+% In statistics, Fisher's method,[1][2] also known as Fisher's combined probability test,
+% is a technique for data fusion or "meta-analysis" (analysis of analyses).
+% It was developed by and named for Ronald Fisher. In its basic form, it is used
+% to combine the results from several independent tests bearing upon the same overall hypothesis (H0).
+if nargin==0
+ disp('usage: p=fishercombine(ps) %add vector with ps of the independent tests');
+end
+chi2=-2*sum(log(ps));
+df=2*length(ps);
+p=1- chi2cdf(chi2,df);
+if nargout==0
+ fprintf('Chi2 = %g df=%d p = %g\n',chi2,df,p);
+else
+ res.p=p;
+ res.chi=chi2;
+ res.df=df;
+end
+end
+
diff --git a/generic_ews.m b/generic_ews.m
new file mode 100644
index 0000000..e2a23d9
--- /dev/null
+++ b/generic_ews.m
@@ -0,0 +1,1008 @@
+%function generic_ews
+% generic_ews is used to estimate statistical moments within rolling
+% windows along a timeserie (based on the R early warning signals toolbox,
+% but simplified)
+%
+% Usage
+%
+% generic_ews(timeseries, option pairs)
+% generic_ews(time,timeseries, option pairs)
+% for example:
+% generic_ews(timeseries, 'winsize', 50, 'detrending', 'no', ...
+% 'bandwidth', [], 'logtransform', false, 'interpolate', false)
+%
+% Arguments
+%
+% timeseries - a numeric vector of the observed univariate timeseries values or a numeric
+% matrix of observed timeseries values. If you use a tabel or dataset
+% you can add time as a column named t or time
+%
+% time - a numeric vector of the time index (numeric)
+%
+% datacolumn - column number with the data (default first column (excluding time))
+%
+% winsize - is the size of the rolling window expressed as percentage of the timeseries length
+% (must be numeric between 0 and 100). Default is 50%.
+%
+% bandwidth - is the bandwidth used for the Gaussian kernel when gaussian filtering is applied.
+% It is expressed as percentage of the timeseries length (must be numeric between 0 and 100)
+% Alternatively it can be given by the optimal bandwidth suggested by Bowman and Azzalini (1997) Default).
+%
+% detrending - the timeseries can be detrended/filtered prior to analysis (can be a cell for each data column).
+% There are four options: 'no'= no detrending, 'gaussian' =
+% gaussian filtering, 'movmean' = moving average (statistics
+% toolbox), 'poly' = polynomical (default degree=4, see
+% polydegree), 'linear' = linear detrending, or 'first-diff' = first-differencing. Default is 'gaussion' detrending.
+%
+% indicators - cell of strings with the names of the indicators can be one of the following strings:
+% 'AR' - autocorrelation of residuals
+% 'acf' - autocorrelation function of residuals
+% 'std' - standard deviation of residuals
+% 'skewness' - skewness of residuals
+% 'cv' - coefficient of variation of original data
+% 'abscorr' - the average of the absolute correlation with all other columns
+% default: {'AR', 'acf', 'std', 'skewness'}
+%
+% figures - cell of strings with figures of the combined figure {'original data', 'residuals';
+% 'indicator 1'...'indicator 4'} 'original data 2' gives the
+% original data of column 2, 'residuals 2' the residuals of
+% col 2, no number = first datacolumn
+%
+% logtransform - if TRUE data are logtransformed prior to analysis as log(X+1). Default is FALSE.
+%
+% interpolate - If TRUE linear interpolation is applied to produce a timeseries of equal length as the original. Default is FALSE (assumes there are no gaps in the timeseries).
+% polydegree - degree of polynomial for detrending (default = 4)
+% arlag - The lag for autoregression. Default = 1
+% silent - If silent=true then the figure is not drawn, default = false
+% ebisuzaki - number of runs for a null model (using the Ebisuzaki method), default = 0
+% violinplot - make violin plot of the taus of the null model, default = true
+% title - use text as title of the figure, default =''
+% corrcolumns - column number with the variable for cross correlation
+% bandwidth_unit - unit of bandwidth (see units of winsize)
+% winsize_unit - units of the winsize - if different from '%' the absolute
+% units of the time series are assumed. If % the real size
+% is calculated from the length of the time series
+% cv - (depreciated) if true use cv instead of standard deviation default = false (use now indicators field)
+%
+%
+%not yet supported in the MATLAB version:
+% AR_n - If TRUE the best fitted AR(n) model is fitted to the data. Default is FALSE.
+%
+% powerspectrum -If TRUE the power spectrum within each rolling window is plotted. Default is FALSE.
+%
+% Options:
+% s=generic_ews('-defaultopts') - create a struct s with the default
+% options (used internally)
+% generic_ews('-f',res) determine the Fisher combined probability plots of
+% the results res. res should be a struct array or cell array of earlier results
+%
+% generic_ews returns a structure with the following fields:
+% timeseries - original time series
+% EWSopt - structure with all used options (including the time index)
+% trend - trend in the time series (depending on the detrending option)
+% colnames - names of columns of the results (indicator names made unique, e.g. if
+% two AR indicators are saved they are named AR and AR_1)
+% indicators - the resulting
+% taus - Kendall tau values for each indicator
+% pvalues - p values of the tau values (only if an ebizuzsaki analysis was
+% done)
+% ebitaus - the tau values for each of the ebizuzsaki iterations
+% done)
+% description - string with main options
+%
+%
+% In addition, generic_ews returns a plots. The plot contains the original data,
+% the detrending/filtering applied and the residuals (if selected), and all the moment statistics.
+% For each statistic trends are estimated by the nonparametric Kendall tau correlation (if statistics toolbox
+% is present).
+% Not supported: The second plot, if asked, quantifies resilience indicators fitting AR(n) selected by the
+% Akaike Information Criterion. The third plot, if asked, is the power spectrum estimated
+% by spec.ar for all frequencies within each rolling window.
+%
+% Author(s)
+%
+% Vasilis Dakos vasilis.dakos@gmail.com
+% MATLAB version by Egbert van Nes
+%
+% References
+%
+% Ives, A. R. (1995). "Measuring resilience in stochastic systems." Ecological Monographs 65: 217-233
+%
+% Dakos, V., et al (2008). "Slowing down as an early warning signal for abrupt climate change." Proceedings of the National Academy of Sciences 105(38): 14308-14312
+%
+% Dakos, V., et al (2012)."Methods for Detecting Early Warnings of Critical Transitions in Time Series Illustrated Using Simulated Ecological Data." PLoS ONE 7(7): e41010. doi:10.1371/journal.pone.0041010
+%
+%
+function result = generic_ews(timeseries, varargin)
+ if nargin == 0
+ error('grind:generic_ews:notenoughargs', 'Not enough input arguments, at least the timeseries is needed (see <a href="matlab:help(''generic_ews'')">help generic_ews</a>)');
+ end
+ if ischar(timeseries)
+ if strcmp(timeseries, '-defaultopts')
+ %the default options
+ result = struct('winsize', 50, ...
+ 'detrending', 'gaussian', ... % = c("no", "gaussian", "linear", "poly", "first-diff", "movmean",
+ 'bandwidth', 10, ...
+ 'time', [], ...
+ 'datacolumn', 1, ... %per indicator or one for all
+ 'corrcolumns', [], .... %empty is all (except datacolumn)
+ 'silent', false, ...
+ 'logtransform', false, ...
+ 'arlag', 1, ...
+ 'indicators', {{'AR', 'acf', 'std', 'skewness'}}, ... %other options cv, abscorr, crosscorr, mean
+ 'figures', {{'original data', 'residuals'; 'indicator 1', 'indicator 2'; 'indicator 3', 'indicator 4'; 'indicator 5', 'indicator 6'}}, ...
+ 'cv', false, ...
+ 'polydegree', 4, ...
+ 'interpolate', false, ...
+ 'AR_n', false, ...
+ 'powerspectrum', false, ...
+ 'bandwidth_unit', '%', ...
+ 'winsize_unit', '%', ...
+ 'ebisuzaki', 0, ... %number of nullmodel (Ebisuzaki) runs (0=none)
+ 'violinplot', true, ...
+ 'title', '');
+ return;
+ elseif strncmpi(timeseries, '-f', 2)
+ %table with fisher combined probabilities
+ res = varargin{1};
+ if ~iscell(res)
+ %if not a cell make a cell for convenience
+ res1 = varargin{1};
+ res = cell(size(res));
+ for i = 1:numel(res1)
+ res{i} = res1(i);
+ end
+ end
+ %ebitaus = zeros(res{1}.EWSopt.ebisuzaki, length(res));
+ titles = cell(length(res), 1);
+ colnames = res{1}.colnames;
+ pvalues = zeros(length(res), length(res{1}.colnames));
+ fisher(1, length(res{1}.colnames)) = struct('p', [], 'chi', [], 'df', []);
+ % fisher = cell(1, length(res{1}.colnames));
+ taus = zeros(length(res), length(res{1}.colnames));
+ silent = res{1}.EWSopt.silent | nargout > 0;
+ for j = 1:length(res{1}.colnames)
+ if ~silent
+ fprintf( '\n======= %s ========\n', res{1}.colnames{j});
+ fprintf( 'Title\tKendall-tau %s\tp-value %s\n', res{1}.colnames{j}, res{1}.colnames{j});
+ end
+ for i = 1:numel(res)
+ %ebitaus(:, i) = res{i}.ebitaus(j, :)';
+ titles{i} = res{i}.EWSopt.title;
+ if isempty(titles{i})
+ titles{i} = sprintf('data set %d', i);
+ end
+ taus(i, j) = res{i}.taus(j);
+ pvalues(i, j) = res{i}.pvalues(j); %(sum(taus(i) < ebitaus(:, i)) + 1) / size(ebitaus, 1);
+ if ~silent
+ fprintf( '%s\t%.3g\t%.3g\n', titles{i}, taus(i, j), pvalues(i, j));
+ end
+ end
+ fisher1 = fishercombine(pvalues(:, j));
+ if ~silent
+ regimes = 'Fisher combined'; %sprintf('%s, ', titles{:});
+ fprintf( '%s\tX2(%d)=%.3g\t%.3g\n', regimes, fisher1.df, fisher1.chi, fisher1.p);
+ end
+ fisher(1, j) = fisher1;
+ end
+ if nargout > 0
+ result.titles = titles;
+ result.colnames = colnames;
+ result.taus = taus;
+ result.pvalues = pvalues;
+ % result.fisherp = fisherp;
+ result.fisher = fisher;
+ end
+ return;
+ end
+ end
+ if ~((isfield(timeseries, 'timeseries') || isfield(timeseries, 'inds')) && isfield(timeseries, 'taus')) % if these fields exist it is a previously saved run
+ [timeseries, EWSopt] = initialize_settings(timeseries, varargin);
+
+ taus = zeros(length(EWSopt.indicators), 1);
+ trend = zeros(size(timeseries));
+ for i = 1:size(timeseries, 2)
+ detrending = getcolumn(EWSopt.detrending, i, true);
+ if strcmpi(detrending, 'first-diff')
+ timeseries(:, i) = [nan; diff(timeseries(:, i))];
+ elseif strcmpi(detrending, 'linear')
+ p = polyfit(EWSopt.time, timeseries(:, i), 1);
+ trend(:, i) = polyval(p, EWSopt.time);
+ elseif strcmpi(detrending, 'poly')
+ p = polyfit(EWSopt.time, timeseries(:, i), EWSopt.polydegree);
+ trend(:, i) = polyval(p, EWSopt.time);
+ elseif strcmpi(detrending, 'movmean')
+ bw = getcolumn(EWSopt.bandwidth, i);
+ if isempty(bw)
+ bw = [4 0];
+ end
+ try
+ trend(:, i) = movmean( timeseries(:, i), bw);
+ catch
+ trend(:, i) = movmean1( timeseries(:, i), bw);
+ end
+ elseif strcmpi(detrending, 'gaussian')
+ bw = getcolumn(EWSopt.bandwidth, i);
+ if strcmp(EWSopt.bandwidth_unit{i}, '%')
+ if ~isempty(bw)
+ EWSopt.absbandwidth = bw / 100 * (max(EWSopt.time) - min(EWSopt.time));
+ else
+ EWSopt.absbandwidth = [];
+ end
+ else
+ if iscell(EWSopt.bandwidth)
+ EWSopt.absbandwidth = EWSopt.bandwidth{i};
+ else
+ EWSopt.absbandwidth=EWSopt.bandwidth;
+ end
+ end
+ trend(:, i) = ksmooth(EWSopt.time, timeseries(:, i), getcolumn(EWSopt.absbandwidth, i), detrending);
+ end
+ end
+
+
+ % ndx = ~isnan(trend(:, 1));
+ if strcmp(EWSopt.winsize_unit, '%')
+ if EWSopt.winsize > 100
+ EWSopt.winsize = 100;
+ end
+ if EWSopt.winsize < 0.1
+ EWSopt.winsize = 0.1;
+ end
+ EWSopt.avwin = EWSopt.winsize / 100 * (max(EWSopt.time) - min(EWSopt.time));
+ else
+ EWSopt.avwin = EWSopt.winsize;
+ end
+
+ indicators = cell(size(EWSopt.indicators));
+ for i1 = 1:length(indicators)
+ indicators{i1} = zeros(size(timeseries, 1), 1) .* NaN;
+ end
+ ts = EWSopt.time;
+ winsize = EWSopt.avwin;
+ tstart = ts - winsize;
+ f = find(tstart >= EWSopt.time(1), 1);
+ if isempty(f)
+ error('not enough data')
+ end
+ tstart = tstart - (tstart(f) - EWSopt.time(1)); %allign with time lags
+ tstart(tstart < EWSopt.time(1)) = NaN;
+ notrend = false;
+ for i1 = 1:length(ts)
+ if ~isnan(tstart(i1)) && ~isnan(trend(i1, 1))
+ ndx = ts >= tstart(i1) & ts < ts(i1) & ~isnan(trend(i1, 1));
+ tim = timeseries(ndx, :);
+ if isempty(trend)
+ tren = zeros(size(tim));
+ else
+ tren = trend(ndx, :);
+ end
+ for j = 1:length(EWSopt.indicators)
+ switch EWSopt.indicators{j}
+ case 'AR'
+ indicators{j}(i1) = autoregression(tim(:, getcolumn(EWSopt.datacolumn, j)) - tren(:, getcolumn(EWSopt.datacolumn, j)), getcolumn(EWSopt.arlag, j));
+ case 'acf'
+ indicators{j}(i1) = acf(tim(:, getcolumn(EWSopt.datacolumn, j)) - tren(:, getcolumn(EWSopt.datacolumn, j)), getcolumn(EWSopt.arlag, j));
+ case 'cv'
+ indicators{j}(i1) = std(tim(:, getcolumn(EWSopt.datacolumn, j))) / mean(tim(:, getcolumn(EWSopt.datacolumn, j)));
+ case 'std'
+ indicators{j}(i1) = std(tim(:, getcolumn(EWSopt.datacolumn, j)) - tren(:, getcolumn(EWSopt.datacolumn, j)));
+ case 'skewness'
+ indicators{j}(i1) = skewness(tim(:, getcolumn(EWSopt.datacolumn, j)) - tren(:, getcolumn(EWSopt.datacolumn, j)));
+ % case 'abscorr'
+ % indicators{j}(i1) = meancorr(tim, getcolumn(EWSopt.datacolumn, j), getcolumn(EWSopt.corrcolumns, j, true), true);
+ % case 'crosscorr'
+ % indicators{j}(i1) = meancorr(tim, getcolumn(EWSopt.datacolumn, j), getcolumn(EWSopt.corrcolumns, j, true), false);
+ % case 'mean'
+ % indicators{j}(i1) = mean(tim(:, getcolumn(EWSopt.datacolumn, j)));
+ case {'abscorr' 'crosscorr', 'mean'}
+ notrend = true;
+ otherwise
+ error('grind:generic_ews:unknown', 'Unknown indicator : %s', EWSopt.indicators{j});
+ end
+ end
+ % fprintf('%d-%d %g %g %g\n',find(ndx,1),find(ndx,1,'last'),ts(find(ndx,1)),ts(find(ndx,1,'last')),indicators{1}(i1))
+ end
+ end
+ if notrend
+ for i1 = 1:length(ts)
+ if ~isnan(tstart(i1))
+ ndx = ts >= tstart(i1) & ts < ts(i1) ;
+ tim = timeseries(ndx, :);
+ for j = 1:length(EWSopt.indicators)
+ switch EWSopt.indicators{j}
+ case 'abscorr'
+ indicators{j}(i1) = meancorr(tim, getcolumn(EWSopt.datacolumn, j), getcolumn(EWSopt.corrcolumns, j, true), true);
+ case 'crosscorr'
+ indicators{j}(i1) = meancorr(tim, getcolumn(EWSopt.datacolumn, j), getcolumn(EWSopt.corrcolumns, j, true), false);
+ case 'mean'
+ indicators{j}(i1) = mean(tim(:, getcolumn(EWSopt.datacolumn, j)));
+ end
+ end
+ % fprintf('%d-%d %g %g %g\n',find(ndx,1),find(ndx,1,'last'),ts(find(ndx,1)),ts(find(ndx,1,'last')),indicators{1}(i1))
+ end
+ end
+ end
+ hasstats = exist('corr', 'file') ~= 0;
+ if hasstats
+ %we cannot use p values as the data are not independent
+ taus = zeros(size(EWSopt.indicators));
+ for i = 1:length(taus)
+ taus(i) = corr(EWSopt.time, indicators{i}, 'type', 'Kendall', 'rows', 'complete');
+ end
+ end
+ ebitaus = [];
+ else %plot previous data
+ hasstats = exist('corr', 'file') ~= 0;
+ res = timeseries;
+ EWSopt = res.EWSopt;
+ EWSopt = mergeoptions(EWSopt, varargin);
+ if ~isfield(EWSopt, 'indicators')
+ EWSopt.indicators = {'AR', 'acf', 'std', 'skewness'};
+ if EWSopt.cv
+ EWSopt.indicators{3} = 'cv';
+ end
+ end
+ % if isa(res.inds, 'table')
+ % inds = table2array(res.inds(:, makeunique(EWSopt.indicators)));
+ % else
+ % inds = res.inds(:, 4:end);
+ % end
+ taus = res.taus;
+ ebitaus = res.ebitaus;
+ if ~isfield(res, 'inds')
+ trend = res.trend;
+ if isempty(trend)
+ trend = zeros(size(res.timeseries));
+ end
+ timeseries = res.timeseries;
+ indicators = num2cell(res.indicators, 1);
+ elseif isa(res.inds, 'table') %old version
+ if ~any(strcmp(res.inds.Properties.VariableNames, 'trend'))
+ trend = zeros(size(res.inds.Variable));
+ else
+ trend = res.inds.trend;
+ end
+ timeseries = res.inds.Variable;
+ inds = table2array(res.inds(:, makeunique(EWSopt.indicators)));
+ indicators = num2cell(inds, 1);
+ end
+ end
+
+
+ pvalues = [];
+
+ if EWSopt.ebisuzaki > 0 && isempty(ebitaus)
+ nulldata = zeros(size(timeseries, 1), size(timeseries, 2), EWSopt.ebisuzaki);
+ for i = 1:size(timeseries, 2)
+ nulldata(:, i, :) = ebisuzaki(timeseries(:, i), EWSopt.ebisuzaki);
+ end
+ %nulldata = ebisuzaki(timeseries(:, EWSopt.datacolumn), EWSopt.ebisuzaki);
+ EWSopt1 = EWSopt;
+ EWSopt1.ebisuzaki = 0;
+ EWSopt1.silent = true;
+ if isempty(ebitaus)
+ ebitaus = zeros(numel(taus), EWSopt.ebisuzaki);
+ for i = 1:EWSopt.ebisuzaki
+ res = generic_ews(nulldata(:, :, i), EWSopt1);
+ ebitaus(:, i) = res.taus;
+ end
+ end
+ end
+ if ~isempty(ebitaus)
+ pvalues = zeros(size(EWSopt.indicators));
+ nebisuzaki = size(ebitaus, 2);
+ for j = 1:length(pvalues)
+ pvalues(j) = (sum(taus(j) < ebitaus(j, :)) + 1) / nebisuzaki;
+ end
+ fprintf('Surrogate data analysis %s\n', EWSopt.title);
+ fprintf('Number of surrogate data sets: %d\n', nebisuzaki);
+ for i = 1:length(pvalues)
+ fprintf('Tau %s = %g, pvalue = %g\n', EWSopt.indicators{i}, taus(i), pvalues(i));
+ end
+ if ~EWSopt.silent && (~isfield(EWSopt, 'violinplot') || EWSopt.violinplot)
+ oldfig = get(0, 'currentfigure');
+ figure;
+ try
+ violinplot(ebitaus', (1:length(pvalues)), 'bandwidth', 0.1, 'categorynames', EWSopt.indicators);
+ catch
+ end
+ adjustplot('-d')
+ hold on
+ plot((1:length(pvalues)), taus, 'ro');
+ ylabel('Kendall \tau')
+ title('Surrogate data')
+ if ~isempty(oldfig)
+ figure(oldfig);
+ end
+ end
+ end
+ if ~EWSopt.silent && ~isempty(EWSopt.figures)
+ % if ~isfield(EWSopt, 'figures')
+ % EWSopt.figures = {'original data', 'residuals'; 'indicator 1', 'indicator 2'; 'indicator 3', 'indicator 4'};
+ % end
+ xlims = [min(EWSopt.time), max(EWSopt.time)];
+ figure;
+ figs = EWSopt.figures'; %strange of MATLAB that this is different
+ for i = 1:numel(EWSopt.figures)
+ h = subplot(size(EWSopt.figures, 1), size(EWSopt.figures, 2), i);
+ if strncmp(figs{i}, 'original data', 13)
+ col = str2double(regexp(figs{i}, '[0-9]*', 'match', 'once'));
+ if isnan(col)
+ col = getcolumn(EWSopt.datacolumn, 1);
+ end
+ plot(EWSopt.time, timeseries(:, col), 'k-', 'Tag', figs{i});
+ xlim(xlims)
+ adjustticklabels(gca, 'Y');
+ if ~isempty(trend)
+ hold on;
+ plot(EWSopt.time, trend(:, col), 'r-', 'Tag', 'trend');
+ hold off;
+ end
+ if ~isempty(EWSopt.title)
+ text(0.1, 0.9, EWSopt.title, 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'toptext')
+ set(gcf, 'name', EWSopt.title)
+ end
+ xlabel('Time');
+ % ylabel('variable and trend');
+ elseif strncmp(figs{i}, 'residuals', 9)
+ col = str2double(regexp(figs{i}, '[0-9]*', 'match', 'once'));
+ if isnan(col)
+ col = getcolumn(EWSopt.datacolumn, 1);
+ end
+ if ~isempty(trend)
+ h = stem(EWSopt.time, timeseries(:, col) - trend(:, col), 'k.', 'Tag', figs{i});
+ adjustticklabels(gca, 'Y');
+ text(0.1, 0.9, 'residuals', 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'toptext')
+ set(h, 'markersize', 1);
+ else
+ plot(EWSopt.time, timeseries(:, col), 'k-', 'Tag', figs{i});
+ text(0.1, 0.9, 'No residuals - no detrending', 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'text1')
+ end
+ xlabel('Time');
+ xlim(xlims);
+ else
+ uinds = makeunique(EWSopt.indicators);
+ iind = find(strcmp(figs{i}, uinds));
+ if isempty(iind)
+ iind = str2double(regexp(figs{i}, '[0-9]*', 'match', 'once'));
+ end
+ if iind <= length(EWSopt.indicators)
+ sindicator = uinds{iind};
+ if any(strcmp(EWSopt.indicators{iind}, {'AR', 'acf'}))
+ sindicator = sprintf('%s(%d)', lower(sindicator), EWSopt.arlag);
+ elseif strcmp(EWSopt.indicators{iind}, 'std')
+ sindicator = 'standard deviation';
+ end
+ plot(EWSopt.time, indicators{iind}, 'k-', 'tag', uinds{iind});
+ adjustticklabels(gca, 'Y');
+ text(0.1, 0.9, sindicator, 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'toptext')
+ if hasstats
+ if ~isempty(pvalues)
+ text(0.1, 0.1, sprintf('{\\tau} = %5.3g (p =%5g)', taus(iind), pvalues(iind)), 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'bottomtext');
+ else
+ text(0.1, 0.1, sprintf('{\\tau} = %5.3g', taus(iind)), 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'bottomtext');
+ end
+ end
+ xlabel('Time');
+ xlim(xlims);
+ else
+ delete(h);
+ end
+ end
+ movegui(gcf, 'center');
+ end
+ ch = findobj(gcf, 'type', 'axes');
+ set(ch, 'TickDir', 'out');
+ if numel(figs) > 1
+ adjustpos([0, 0.07]);
+ removeaxis('x');
+ end
+ xlim(xlims);
+ end
+ if nargout > 0
+ % indic = makeunique(EWSopt.indicators);
+ % if ~isoctave && verLessThan('matlab', '8.4.0')
+ % properties = struct('VarNames', indic , 'Description', sprintf('data generated with generic_ews\nwinsize = %g, detrending = %s, bandwidth = %g, logtransform = %d, arlag = %d, interpolate = %d', EWSopt.winsize, getcolumn(EWSopt.detrending, 1), getcolumn(EWSopt.bandwidth, 1), EWSopt.logtransform, EWSopt.arlag, EWSopt.interpolate));
+ % inds = struct('time', EWSopt.time, 'timeseries', timeseries, 'trend', trend, 'indicators', [indicators{:}], 'Properties', properties);
+ % else
+ % % [~,ndx]=unique(EWSopt.indicators);
+ % if ~isempty(trend)
+ % ndx = ~isnan(trend(:, 1));
+ % for i = 1:length(indicators)
+ % inds = zeros(size(ndx)) + NaN;
+ % if numel(indicators{i}) == numel(ndx)
+ % inds = indicators{i};
+ % else
+ % inds(ndx) = indicators{i};
+ % end
+ % indicators{i} = inds;
+ % end
+ % inds = table(EWSopt.time, timeseries, trend, indicators{:});
+ % inds.Properties.VariableNames = [{'Time', 'Variable', 'trend'}, indic];
+ % else
+ % inds = table(EWSopt.time, timeseries, indicators{:});
+ % inds.Properties.VariableNames = [{'Time', 'Variable'}, indic];
+ % end
+ %
+ % inds.Properties.Description = sprintf('data generated with generic_ews\nwinsize = %g, detrending = %s, bandwidth = %g, logtransform = %d, arlag = %d, interpolate = %d', EWSopt.winsize, getcolumn(EWSopt.detrending, 1), getcolumn(EWSopt.bandwidth, 1), EWSopt.logtransform, EWSopt.arlag, EWSopt.interpolate);
+ % end
+ if isempty(getcolumn(EWSopt.bandwidth, 1))
+ % optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
+ n = size(timeseries, 1);
+ hx = median(abs(EWSopt.time - median(EWSopt.time))) / 0.6745 * (4 / 3 / n)^0.2;
+ hy = median(abs(timeseries(:, getcolumn(EWSopt.datacolumn, 1)) - median(timeseries(:, getcolumn(EWSopt.datacolumn, 1))))) / 0.6745 * (4 / 3 / n)^0.2;
+ EWSopt.bandwidth = sqrt(hy * hx);
+ end
+ if nargout == 1
+ descr = sprintf('data generated at %s with generic_ews\nwinsize = %g, detrending = %s, bandwidth = %g, logtransform = %d, arlag = %d, interpolate = %d', ...
+ datestr(now()), EWSopt.winsize, getcolumn(EWSopt.detrending, 1), getcolumn(EWSopt.bandwidth, 1), EWSopt.logtransform, EWSopt.arlag, EWSopt.interpolate);
+ result = struct('timeseries', timeseries, ...
+ 'EWSopt', EWSopt, ...
+ 'trend', trend, ...
+ 'colnames', {makeunique(EWSopt.indicators)} , ...
+ 'indicators', [indicators{:}], ...
+ 'taus', taus, ...
+ 'pvalues', pvalues, ...
+ 'ebitaus', ebitaus, ...
+ 'description', descr);
+ end
+ end
+
+ %end
+end
+function [timeseries, EWSopt] = initialize_settings(timeseries, args)
+ if isa(timeseries, 'table')
+ %the time series may be a table , it can then have a column with
+ %time, the name of that column should be t or Time or time
+ ndx = ~(strcmp('t', timeseries.Properties.VariableNames) | strcmpi('time', timeseries.Properties.VariableNames));
+ if any(ndx)
+ % the time is added to the options
+ args = [args, 'time', timeseries.t];
+ timeseries = timeseries{:, ndx};
+ else
+ timeseries = timeseries{:, :};
+ end
+ elseif isa(timeseries, 'dataset')
+ %the time series may be a dataset, it can then have a column with
+ %time, the name of that column should be t or Time or time
+ ndx = ~(strcmp('t', timeseries.Properties.VarNames) | strcmpi('time', timeseries.Properties.VarNames));
+ if any(ndx)
+ % the time is added to the options
+ args = [args, 'time', timeseries.t];
+ timeseries = double(timeseries(:, ndx));
+ else
+ timeseries = double(timeseries(:, :));
+ end
+ else
+ timeseries = double(timeseries);
+ end
+
+% if ~exist('i_use', 'file')
+% addpath([grindpath, filesep, 'sys2']);
+% end
+ if isempty(timeseries)
+ error('grind:generic_ews:empty', 'Empty time series');
+ end
+ EWSopt = generic_ews('-defaultopts');
+ EWSopt = mergeoptions(EWSopt, args);
+ if ischar(EWSopt.bandwidth_unit)
+ EWSopt.bandwidth_unit={EWSopt.bandwidth_unit};
+ end
+ if numel(EWSopt.bandwidth_unit)==1&&size(timeseries,2)>1
+ EWSopt.bandwidth_unit=repmat(EWSopt.bandwidth_unit,1,size(timeseries,2));
+ end
+
+ if EWSopt.powerspectrum || EWSopt.AR_n
+ error('Option not yet supported');
+ end
+ if isempty(EWSopt.time)
+ EWSopt.time = transpose(linspace(1, size(timeseries, 1), size(timeseries, 1)));
+ end
+
+ if EWSopt.interpolate
+ ts = transpose(linspace(EWSopt.time(1), EWSopt.time(end), size(timeseries, 1)));
+ timeseries = interp1(EWSopt.time, timeseries, ts);
+ EWSopt.time = ts;
+ end
+ for i = 1:size(timeseries, 2)
+ if getcolumn(EWSopt.logtransform, i)
+ timeseries(:, i) = log(timeseries(:, i) + 1);
+ end
+ end
+
+ EWSopt.datalength = max(EWSopt.time) - min(EWSopt.time);
+ if EWSopt.cv
+ %depreciated
+ f = strcmp(EWSopt.indicators, 'std');
+ EWSopt.indicators(f) = {'cv'};
+ end
+end
+function EWSopt = mergeoptions(EWSopt, arg)
+ f = fieldnames(EWSopt);
+ if ~isempty(arg) && isstruct(arg{1})
+ %all options can be entered as struct
+ if exist('mergestructs', 'file') == 0
+ %initgrind;
+ end
+ %replace common fields in EWSopt
+ EWSopt = mergestructs(EWSopt, arg{1});
+ else
+ for i = 1:2:length(arg)
+ if ~any(strcmp(f, arg{i}))
+ s = sprintf('"%s", ', f{:});
+ error('Unknown option for generic_ews: %s\nValid options are: %s', arg{i}, s(1:end - 1));
+ end
+ EWSopt.(arg{i}) = arg{i + 1};
+ end
+ end
+end
+
+
+%merge two structs
+%astruct=mergestructs(astruct,newstruct)
+%
+%fields of a astruct will be replaced by newstruct, new fields will be
+%added if necessary
+%nested for substructures
+function fullstruct = mergestructs(fullstruct, substruct)
+ if isempty(fullstruct)
+ fullstruct = substruct;
+ elseif ~isempty(substruct) && numel(substruct) == numel(fullstruct)
+ fnew = fieldnames(substruct);
+ fold = fieldnames(fullstruct);
+ for i = 1:length(fnew)
+ for j = 1:numel(substruct)
+ if isstruct(substruct(j).(fnew{i})) && any(strcmp(fold, fnew{i})) && isstruct(fullstruct(j).(fnew{i}))
+ fullstruct(j).(fnew{i}) = mergestructs(fullstruct(j).(fnew{i}), substruct(j).(fnew{i}));
+ else
+ fullstruct(j).(fnew{i}) = substruct(j).(fnew{i});
+ end
+
+ end
+
+ end
+
+ else
+ error('grind:mergestructs', 'Cannot merge strucs of different sizes')
+ end
+end
+
+
+function res = getcolumn(col, i, multcell)
+ if nargin < 3
+ multcell = false;
+ end
+ if multcell && ~iscell(col)
+ col = {col};
+ end
+ if isempty(col)
+ res = col;
+ elseif iscell(col)
+ if numel(col) == 1
+ res = col{1};
+ else
+ res = col{i};
+ end
+ else
+ if numel(col) == 1
+ res = col(1);
+ else
+ res = col(i);
+ end
+ end
+end
+
+function res = meancorr(x, datacolumm, corrcolumns, doabs)
+ if size(x, 2) == 1
+ res = NaN;
+ return;
+ end
+ res = corr(x);
+ res = res(datacolumm, :);
+ if isempty(corrcolumns)
+ res(:, datacolumm) = [];
+ else
+ res = res(:, corrcolumns);
+ end
+ if doabs
+ res = mean(abs(res));
+ else
+ res = mean(res);
+ end
+end
+function indic = makeunique(indic)
+ uindic = unique(indic);
+ if length(uindic) < length(indic)
+ for i = 1:length(uindic)
+ f = find(strcmp(uindic{i}, indic));
+ for j = 2:length(f)
+ indic{f(j)} = sprintf('%s_%d', indic{f(j)}, j - 1);
+ end
+ end
+ end
+end
+function b = linregres(x, y) %#ok<DEFNU>
+ %b=(sum(x.*y)-1/length(x)*sum(x)*sum(y))/(sum(x.^2)-1/length(x)*sum(x).^2);
+ b = sum(x .* y) ./ (sum(x.^2)); %without intercept
+ %intercept
+ %a=mean(y)-b*mean(x);
+end
+function res = skewness(x, flag)
+ if nargin < 2
+ flag = 0;
+ end
+ x = x - mean(x);
+ s2 = mean(x.^2); % this is the biased variance estimator
+ m3 = mean(x.^3);
+ res = m3 ./ s2.^(1.5);
+ if flag == 0
+ n = length(x);
+ if n > 3
+ n(n < 3) = NaN; % bias correction is not defined for n < 3.
+ res = res .* sqrt((n - 1) ./ n) .* n ./ (n - 2);
+ end
+ end
+end
+function res = acf(X, alag)
+ if nargin < 2
+ alag = 1;
+ end
+ %autocorrelation function like in R
+ n = length(X);
+ s = var(X);
+ mu = mean(X);
+ Xt = X(1:end - alag);
+ Xtk = X(1 + alag:end);
+ res = 1 / (n - 1) / s * sum((Xt - mu) .* (Xtk - mu));
+end
+function res = ksmooth(x, y, bandwidth, method)
+ if nargin < 3
+ bandwidth = [];
+ end
+ if nargin < 4
+ method = 'g';
+ end
+ if isempty(x)
+ res = x;
+ return
+ end
+ n = length(x);
+ if isempty(bandwidth)
+ % optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
+ hx = median(abs(x - median(x))) / 0.6745 * (4 / 3 / n)^0.2;
+ hy = median(abs(y - median(y))) / 0.6745 * (4 / 3 / n)^0.2;
+ h = sqrt(hy * hx);
+ fprintf('bandwidth smoothing set to %g\n', h);
+ if h < sqrt(eps) * n
+ error('GRIND:outf:LittleVar', 'There is not enough variation in the data. Regression is meaningless.')
+ end
+ else
+ h = bandwidth;
+ end
+ switch lower(method(1))
+ case 'g' %"Normal" = Gaussian kernel function
+ h = h * 0.36055512754640; %variance of 0.13 to get quartiles at + / - 0.25 * h (comparable with ksmooth (R))
+ res = ones(size(y, 1), 1);
+ for k = 1:n
+ xx = abs(x(k) - x) / h;
+ z = exp( - xx(xx < 4).^2 / 2); %x < 4 is more efficient for long time series (negligible effect)
+ res(k) = sum(z .* y(xx < 4)) / sum(z);
+ end
+ %remove tails
+ x1 = (x(end) - x) > h / 2 & (x - x(1)) > h / 2;
+ res(~x1) = NaN;
+ case 'b' %"box" = moving average
+ d = h / 2; % 0.25 quartiles
+ res = cell(1, n);
+ for k = 1:n
+ xx = abs(x(k) - x) / h;
+ z = xx < d;
+ res(k) = sum(z .* y) / sum(z);
+ end
+ end
+end
+
+function res = movmean1( timeseries, bw)
+ %simple remake of movmean (for older MAT
+ if length(bw) == 1
+ if mod(bw, 2) == 0
+ bw = [bw / 2 bw / 2 - 1];
+ else
+ bw = [(bw - 1) / 2 (bw - 1) / 2];
+ end
+ end
+ starti = [1:length(timeseries)] - bw(1);
+ endi = [1:length(timeseries)] + bw(2);
+ starti(starti < 1) = 1;
+ endi(endi > length(timeseries)) = length(timeseries);
+ res = zeros(size(timeseries));
+ for i1 = 1:length(timeseries)
+ res(i1) = mean(timeseries(starti(i1):endi(i1)));
+ end
+end
+
+function removeaxis(orientation)
+ hs = getsubplotgrid(gcf);
+ if nargin == 0
+ orientation = 'b';
+ end
+ [rows, cols] = size(hs);
+ if strncmpi(orientation, 'b', 1) || strncmpi(orientation, 'x', 1)
+ for i = 1:rows - 1
+ for j = 1:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ set(get(hs(i, j), 'xlabel'), 'string', '')
+ set(hs(i, j), 'XTickLabel', []);
+ end
+ end
+ end
+ for i = 2:rows
+ for j = 1:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ set(get(hs(i, j), 'title'), 'string', '')
+ end
+ end
+ end
+ end
+ if strncmpi(orientation, 'b', 1) || strncmpi(orientation, 'y', 1)
+ for i = 1:rows
+ for j = 2:cols
+ if ishandle(hs(i, j)) && (hs(i, j) ~= 0)
+ set(get(hs(i, j), 'ylabel'), 'string', '')
+ set(hs(i, j), 'YTickLabel', []);
+ end
+ end
+ end
+ end
+ set(hs(hs ~= 0), 'fontsize', 12)
+ aa = get(get(gcf, 'children'), 'xlabel');
+ if ~iscell(aa)
+ aa = {aa};
+ end
+ for i = 1:length(aa)
+ set(aa{i}, 'fontsize', 12);
+ end
+ aa = get(get(gcf, 'children'), 'ylabel');
+ if ~iscell(aa)
+ aa = {aa};
+ end
+ for i = 1:length(aa)
+ set(aa{i}, 'fontsize', 12);
+ end
+end
+function adjustpos(spacing)
+ set(gcf, 'PaperPositionMode', 'auto')
+ set(gcf, 'units', 'normalized');
+ figpos = get(gcf, 'position');
+ vert = 0.5469; %normal size of figure
+ hor = 0.41;
+ hs = getsubplotgrid(gcf);
+ if nargin == 0
+ spacing = 0.02;
+ end
+ if length(spacing) == 1
+ spacing = spacing + zeros(2, 1);
+ end
+ [rows, cols] = size(hs);
+ if cols < 4 && rows <= 4
+ set(gcf, 'position', [figpos(1:2) hor * cols / 4 vert * rows / 4]);
+ else
+ rat = rows / cols;
+ if rat < 1
+ set(gcf, 'position', [figpos(1:2) hor vert * rat]);
+ else
+ set(gcf, 'position', [figpos(1:2) hor / rat vert]);
+ end
+ end
+ for i = 1:rows
+ for j = 1:cols
+ h = findobj(gcf, 'tag', sprintf('subplot(%d,%d)', i, j));
+ if any(ishandle(h)) && any(h ~= 0)
+ set(h, 'position', subplotpos(rows, cols, j, i, spacing));
+ end
+ end
+ end
+end
+
+function pos = subplotpos(rows, cols, x1, y1, spacing, colwidth)
+ %tight position for subplots, you can only set the spacing between the
+ %figures (simplified from subaxis)
+ % s=subplot('position',subplotpos(rows,cols,x1,y1),'tag',sprintf('subplot(%d,%d)',x1,y1));
+ %
+ if nargin < 5
+ spacing = [0.02 0.02];
+ end
+ if nargin<6
+ colwidth=1;
+ end
+ if length(spacing) == 1
+ spacing = spacing + zeros(2, 1);
+ end
+
+ Args = struct('Holdaxis', 0, ...
+ 'SpacingVertical', spacing(1), 'SpacingHorizontal', spacing(2), ...
+ 'MarginLeft', .2, 'MarginRight', .1, 'MarginTop', 0.05, 'MarginBottom', .1, ...
+ 'rows', rows, 'cols', cols);
+
+
+ cellwidth = ((1 - Args.MarginLeft - Args.MarginRight) - (Args.cols - 1) * Args.SpacingHorizontal) / Args.cols;
+ if colwidth>1
+ cellwidth=cellwidth*colwidth+(colwidth-1)*Args.SpacingHorizontal;
+ end
+ cellheight = ((1 - Args.MarginTop - Args.MarginBottom) - (Args.rows - 1) * Args.SpacingVertical) / Args.rows;
+ xpos1 = Args.MarginLeft + cellwidth * (x1 - 1) + Args.SpacingHorizontal * (x1 - 1);
+ xpos2 = Args.MarginLeft + cellwidth * x1 + Args.SpacingHorizontal * (x1 - 1);
+ ypos1 = Args.MarginTop + cellheight * (y1 - 1) + Args.SpacingVertical * (y1 - 1);
+ ypos2 = Args.MarginTop + cellheight * y1 + Args.SpacingVertical * (y1 - 1);
+
+ pos = [xpos1 1 - ypos2 xpos2 - xpos1 ypos2 - ypos1];
+end
+
+function hs = getsubplotgrid(h)
+ ch = get(h, 'children');
+ tags = get(ch, 'tag');
+ ch = ch(~strcmp(tags, 'legend'));
+ types = get(ch, 'type');
+ poss = get(ch, 'position');
+ if length(ch) == 1
+ poss = {poss};
+ types = {types};
+ end
+ ipos = zeros(length(ch), 4);
+ for i = length(ch): - 1:1
+ if iscell(types) && strcmp(types{i}, 'axes')
+ ipos(i, :) = poss{i};
+ elseif ischar(types) && strcmp(types, 'axes')
+ ipos(i, :) = poss{i};
+ else
+ ipos(i, :) = [];
+ end
+ end
+ colpos = sort(unique(sort(ipos(:, 1))), 'ascend');
+ rowpos = sort(unique(sort(ipos(:, 2))), 'descend');
+ hs = zeros(length(rowpos), length(colpos));
+ for i = 1:length(ch)
+ if strcmp(types{i}, 'axes')
+ arow = find(rowpos == ipos(i, 2));
+ acol = find(colpos == ipos(i, 1));
+ hs(arow, acol) = ch(i);
+ set(ch(i), 'tag', sprintf('subplot(%d,%d)', arow, acol));
+ end
+ end
+end
+
+%autoregression ols (without intercept and with subtraction the mean)
+%
+function rc = autoregression(Ys, alag)
+ if nargin < 2
+ alag = 1;
+ end
+
+ rc = zeros(length(alag), size(Ys, 2));
+ Ys = Ys - mean(Ys(2:end)); %subtracting the mean
+ % (important if you do regression without intercept!)
+ for j = 1:length(alag)
+ for i = 1:size(Ys, 2)
+ res = linregres(Ys(1:end - alag, i), Ys(1 + alag:end, i));
+ rc(j, i) = res(1);
+ end
+
+ end
+end
+
+
+
+function adjustticklabels(hax, orient)
+ % axis(hax, 'tight')
+ % newtick = get(hax, [orient 'tick']);
+ % tickdiff = (newtick(2) - newtick(1));
+ % newlim = [newtick(1) - tickdiff newtick(end) + tickdiff];
+ % axis(hax, 'manual')
+ % set(hax, [orient 'lim'], newlim);
+ % set(hax, [orient 'tick'], [newtick(1) - tickdiff newtick]);
+end
+
+
+
+
diff --git a/generic_ews_sens.m b/generic_ews_sens.m
new file mode 100644
index 0000000..ce9d45f
--- /dev/null
+++ b/generic_ews_sens.m
@@ -0,0 +1,454 @@
+%function generic_ews_sens
+% make heatmaps of the sensitivity of the tau values of generic_ews of the
+% window size and bandwithd
+%
+% Usage
+%
+% generic_ews_sens(timeseries, option pairs)
+%
+% all options pairs: (with defaults)
+% res= generic_ews_sens(timeseries, 'winsize', (5:5:90), 'detrending', 'no', ...
+% 'bandwidth', (1:1:20), 'logtransform', false, 'interpolate', false)
+% To draw plots:
+% generic_ews_sens('-p',res) plot all results.
+% generic_ews_sens('-p',res,'indicators',{'AR'},'units',unit) - some options
+% can be changed for plotting:
+% - 'indicators' the indicators to be plotted
+% - 'defaultpos' position of the default settings
+% - 'units' the units for the bandwidth and window size in the plot
+% - 'title' title of the figure, default the indicator
+% other changed settings are ignored
+%
+% generic_ews_sens('-f',res) determine the Fisher combined probability plots of
+% the results res. res should be a struct array or cell array of earlier results
+% s=generic_ews_sens('-defaultopts') - create a struct s with the default
+% settings (used internally)
+%
+%
+% Arguments (only those different from generic_ews)
+%
+% winsize - an array with the window sizes to be tested (in % of the
+% size of the time series)
+%
+% bandwidth - an array with the bandwidths to be tested (in % of the
+% size of the time series)
+% defaultpos - position of a star with the default settings (bandwidth,winsize) (optional)
+% bandwidth_unit - unit of bandwidth (see units of winsize)
+% winsize_unit - units of the winsize - if different from '%' the absolute
+% units of the time series are assumed. If % the real size
+% is calculated from the length of the time series
+% distribution - way of selecting bandwidth/winsize (regular/uniform)
+% (default regular)
+% ndraws - number of draws from distribution (ignored when regular)
+% (default = 1000)
+%
+%
+% generic_ews_sens returns a struct with the fields:
+%
+% bandwidths: matrix with all bandwidth tested
+% winsizes: matrix with the window sizes
+% taus: 3D matix with the taus (taus(:,:,1)=AR1,
+% taus(:,:,2)=acf, taus(:,:,3)=standard dev/CV, taus(:,:,4)=skewness)
+% EWSopt the options used
+% inds the timeseries and the lastly calculated indicators (see generic_ews)
+% colnames the names of the columns (see generic_ews)
+%
+%
+%
+% Author(s)
+%
+% Vasilis Dakos vasilis.dakos@gmail.com
+% MATLAB version by Egbert van Nes
+%
+% References
+%
+% Dakos, V., et al (2008). "Slowing down as an early warning signal for abrupt climate change." Proceedings of the National Academy of Sciences 105(38): 14308-14312
+%
+% Dakos, V., et al (2012)."Methods for Detecting Early Warnings of Critical Transitions in Time Series Illustrated Using Simulated Ecological Data." PLoS ONE 7(7): e41010. doi:10.1371/journal.pone.0041010
+%
+%
+function [res] = generic_ews_sens(timeseries, varargin)
+ if isfield(timeseries, 'timeseries')
+ varargin = [{timeseries} varargin];
+ timeseries = '-p';
+ end
+ if nargin > 0 && ischar(timeseries)
+ if strncmp(timeseries, '-p', 2)
+ %plot the figures
+ if numel(varargin{1}) > 1
+ res = varargin{1};
+ for i = 1: numel(varargin{1})
+ generic_ews_sens('-p', res(i), varargin{2:end});
+ end
+ return;
+ end
+ EWSopt = varargin{1}.EWSopt;
+ if length(varargin) > 1
+ EWSopt = mergeoptions(EWSopt, varargin(2:end));
+ end
+ results = varargin{1};
+ %bandwidths = varargin{1}.bandwidths;
+ %winsizes = varargin{1}.winsizes;
+ taus = varargin{1}.taus;
+ ebitaus = varargin{1}.ebitaus;
+ colnames = varargin{1}.colnames;
+ %EWSopt.indicators = varargin{1}.colnames;
+ if ischar( EWSopt.indicators)
+ EWSopt.indicators = {EWSopt.indicators};
+ end
+ figs = EWSopt.indicators;
+ if ischar(figs)
+ figs = {figs};
+ elseif isnumeric(figs)
+ figs = EWSopt.indicators(figs);
+ end
+ if ischar(EWSopt.bandwidth_unit)
+ EWSopt.bandwidth_unit = {EWSopt.bandwidth_unit};
+ end
+ if numel(EWSopt.bandwidth_unit) == 1 && size(timeseries, 2) > 1
+ EWSopt.bandwidth_unit = repmat(EWSopt.bandwidth_unit, 1, size(timeseries, 2));
+ end
+ % bw_unit = EWSopt.bandwidth_unit;
+ % win_unit = EWSopt.winsize_unit;
+ % if isfield(EWSopt, 'units')
+ % units = EWSopt.units;
+ % else
+ % units = {'%', '%'};
+ % end
+ % if ischar(units)
+ % units = {units, units};
+ % end
+ % if length(units) == 1
+ % units = [units units];
+ % end
+ % if any(strcmp(bw_unit{1}, '%'))
+ % bandwidths = bandwidths .* EWSopt.datalength ./ 100;
+ % end
+ % if strcmp(win_unit, '%')
+ % winsizes = winsizes .* EWSopt.datalength ./ 100;
+ % end
+ if isfield(EWSopt, 'defaultpos')
+ defaultpos = EWSopt.defaultpos;
+ if iscell(defaultpos)
+ defaultpos = [defaultpos{:}];
+ end
+ else
+ defaultpos = [];
+ end
+ titles = colnames;
+ for i = 1:numel(figs)
+ fignr = find(strcmp(colnames, figs{i}));
+ if isempty(fignr)
+ error('grind:generic_ews_sens:unknowindicator', 'indicator %s unknown', figs{i});
+ end
+ if ~isempty(EWSopt.title)
+ atitle = EWSopt.title;
+ else
+ atitle = titles{fignr};
+ end
+ plot_sensitiv(results, taus(:, :, fignr), fignr, [], [], defaultpos, atitle, EWSopt);
+ end
+ if ~isempty(ebitaus)
+ if ~isfield(EWSopt,'distribution')||strcmp(EWSopt.distribution, 'regular')
+ titles = colnames;
+ levels = [0:0.05:0.3];
+ for i = 1:numel(figs)
+ fignr = find(strcmp(colnames, figs{i}));
+ if ~isempty(EWSopt.title)
+ atitle = ['p-value''s ' EWSopt.title ];
+ else
+ atitle = ['p-value''s ' titles{fignr}];
+ end
+ pvalues = getpvalues(taus(:, :, fignr), ebitaus(:, :, fignr, :));
+ plot_sensitiv(results, pvalues, fignr, [0 0.2], levels, defaultpos, atitle, EWSopt);
+ end
+ else
+ for i = 1:numel(figs)
+ fignr = find(strcmp(colnames, figs{i}));
+ if ~isempty(EWSopt.title)
+ atitle = ['p-value''s ' EWSopt.title ];
+ else
+ atitle = ['p-value''s ' titles{fignr}];
+ end
+ pvalues = getpvalues(taus(:, :, fignr), ebitaus(:, :, fignr, :));
+ figure;
+ violinplot(pvalues, ones(size(pvalues)));
+ title(atitle)
+ end
+ end
+ end
+ return;
+ elseif strcmp(timeseries, '-setcolor')
+ c1 = findobj(gcf, 'tag', 'perc_fill');
+ if ~isempty(c1)
+ %something wrong with this code... (it works only once)
+ colors = get(c1, 'FaceColor');
+ colors = vertcat(colors{:});
+ [~, i] = max(var(colors, 1));
+ rangecolors = 1 - (colors - max(colors(:, i))) ./ (min(colors(:, i) - max(colors(:, i))));
+ mincolor = varargin{1};
+ if ischar(mincolor)
+ colorcodes = {[1 1 0], 'y'; ...
+ [1 0 1], 'm'; ...
+ [0 1 1], 'c'; ...
+ [1 0 0], 'r'; ...
+ [0 1 0], 'g'; ...
+ [0 0 1], 'b'; ...
+ [1 1 1], 'w'; ...
+ [0 0 0], 'k'};
+ ndx = strcmp(mincolor, colorcodes(:, 2));
+ mincolor = colorcodes{ndx,1};
+ end
+ maxcolor = ones(1, 3);
+ for i = 1:length(c1)
+ set(c1(i), 'facecolor', (mincolor + (maxcolor - mincolor) * (rangecolors(i))));
+ end
+ else
+ error('generic_ews_sens:notfound', 'series not found');
+ end
+ return;
+ elseif strncmpi(timeseries, '-f', 2) %fisher combine previous results
+ inputs1 = varargin{1};
+ inputs1 = generic_ews_sens('-add_pvalues', inputs1);
+ EWSopt = inputs1(1).EWSopt;
+ if length(varargin) > 1
+ EWSopt = mergeoptions(EWSopt, varargin(2:end));
+ end
+ ind_ndx = find(ismember(inputs1(1).colnames, EWSopt.indicators));
+ if isempty(ind_ndx)
+ error('generic_ews_sens:fisher', 'None of the indicators were found');
+ end
+ res1.pvalues = {inputs1(:).pvalues};
+ res1.fisher_pvalues = ones(size(inputs1(1).taus, 1), size(inputs1(1).taus, 2), length(ind_ndx));
+ res1.fisher_chi = ones(size(res1.fisher_pvalues));
+ titles = inputs1(1).colnames;
+ levels = [0:0.05:0.3];
+ if iscell(EWSopt.defaultpos)
+ defaultpos = [EWSopt.defaultpos{:}];
+ else
+ defaultpos = EWSopt.defaultpos;
+ end
+ for k = 1:length(ind_ndx) %loop the different indicators
+ for f1 = 1:size(inputs1(1).taus, 1)
+ for f2 = 1:size(inputs1(1).taus, 2)
+ ps = zeros(size(inputs1));
+ for i = 1:numel(inputs1)
+ ps(i) = inputs1(i).pvalues(f1, f2, ind_ndx(k));
+ end
+ fisher1 = fishercombine(ps);
+ res1.fisher_pvalues(f1, f2, k) = fisher1.p;
+ res1.fisher_chi(f1, f2, k) = fisher1.chi;
+ end
+ end
+ res1.fisher_df = fisher1.df;
+ end
+ res1.hfigs = [];
+ for k = 1:length(ind_ndx)
+ if ~isempty(EWSopt.title)
+ atitle = ['Fisher p ' EWSopt.title ];
+ else
+ atitle = ['Fisher p ' titles{ind_ndx(k)}];
+ end
+ res1.hfigs(k) = plot_sensitiv(inputs1(1), res1.fisher_pvalues(:, :, k), ...
+ k, [0 0.3], levels, defaultpos, atitle, EWSopt);
+ end
+ for i = 1:length(res1.pvalues)
+ res1.pvalues{i} = res1.pvalues{i}(:, :, ind_ndx);
+ end
+ if nargout > 0
+ res = res1;
+ end
+ return;
+ elseif strcmp(timeseries, '-add_pvalues') %get the p-values
+ res = varargin{1};
+ res(1).pvalues = [];
+ if numel(res) > 1
+ for i = numel(res): -1:1
+ res(i) = generic_ews_sens('-add_pvalues', res(i));
+ end
+ else
+ res.pvalues = ones(size(res.taus));
+ for k = 1:size(res(1).taus, 3)
+ res.pvalues(:, :, k) = getpvalues(res.taus(:, :, k), res.ebitaus(:, :, k, :));
+ end
+ end
+ return
+ elseif strcmp(timeseries, '-defaultopts')
+ res = generic_ews('-defaultopts');
+ %default winsizes and bandwidths
+ res.winsize = (5:5:90);
+ res.bandwidth = (1:1:20);
+ res.defaultpos = {[]};
+ res.distribution = 'regular';
+ res.ndraws = 1000;
+ return;
+ end
+ end
+ EWSopt = generic_ews_sens('-defaultopts');
+ EWSopt = mergeoptions(EWSopt, varargin);
+ sens_opts = EWSopt;
+
+ %override the value of silent to suppress too many figures
+ EWSopt.silent = true;
+% if exist('i_waitbar', 'file') == 0
+% initgrind;
+% end
+ if strcmp(EWSopt.distribution, 'regular')
+ inds = [];
+ taus = zeros(numel(EWSopt.winsize), numel(EWSopt.bandwidth), numel(EWSopt.indicators));
+ ebitaus = zeros(numel(EWSopt.winsize), numel(EWSopt.bandwidth), numel(EWSopt.indicators), EWSopt.ebisuzaki);
+ [bandwidths, winsizes] = meshgrid(EWSopt.bandwidth, EWSopt.winsize);
+ % wb = i_parfor_waitbar(0, size(winsizes, 1), 'Calculating sensitivity');
+ parfor i = 1:size(winsizes, 1)
+ taus_loc = taus(i, :, :);
+ ebitaus_loc = ebitaus(i, :, :, :);
+ for j = 1:size(winsizes, 2)
+ EWSopt1 = EWSopt;
+ EWSopt1.winsize = winsizes(i, j);
+ EWSopt1.bandwidth = bandwidths(i, j);
+ res2 = generic_ews(timeseries, EWSopt1);
+ taus_loc(1, j, :) = res2.taus;
+ if ~isempty(res2.ebitaus)
+ ebitaus_loc(1, j, :, :) = res2.ebitaus;
+ end
+ end
+ % i_parfor_waitbar(wb, i);
+ taus(i, :, :) = taus_loc;
+ ebitaus(i, :, :, :) = ebitaus_loc;
+ end
+ elseif strcmp(EWSopt.distribution, 'uniform')
+ taus = zeros(1, EWSopt.ndraws, numel(EWSopt.indicators));
+ inds = zeros(EWSopt.ndraws, size(timeseries, 1), numel(EWSopt.indicators));
+ ebitaus = zeros(1, EWSopt.ndraws, numel(EWSopt.indicators), EWSopt.ebisuzaki);
+ bandwidths = min(EWSopt.bandwidth) + rand(EWSopt.ndraws, 1) * (max(EWSopt.bandwidth) - min(EWSopt.bandwidth));
+ winsizes = min(EWSopt.winsize) + rand(EWSopt.ndraws, 1) * (max(EWSopt.winsize) - min(EWSopt.winsize));
+ %wb = i_parfor_waitbar(0, EWSopt.ndraws, 'Calculating sensitivity');
+ for i = 1:EWSopt.ndraws
+ EWSopt1 = EWSopt;
+ EWSopt1.winsize = winsizes(i);
+ EWSopt1.bandwidth = bandwidths(i);
+ res2 = generic_ews(timeseries, EWSopt1);
+ taus(1, i, :) = res2.taus;
+ if ~isempty(res2.ebitaus)
+ ebitaus(1, i, :, :) = res2.ebitaus;
+ end
+ inds(i, :, :) = res2.indicators;
+ %i_parfor_waitbar(wb, i);
+ end
+ end
+ EWSopt1 = EWSopt;
+ EWSopt1.winsize = winsizes(1);
+ EWSopt1.bandwidth = bandwidths(1);
+ EWSopt1.ebisuzaki = 0;
+
+ inputs1 = generic_ews(timeseries, EWSopt1);
+ EWSopt1.ebisuzaki = EWSopt.ebisuzaki;
+ %i_parfor_waitbar([]);
+ %save the basic data per run
+ if ~isempty(inds)
+ inputs1.indicators = inds;
+ end
+ inputs1.bandwidths = bandwidths;
+ inputs1.winsizes = winsizes;
+ inputs1.taus = taus;
+ if EWSopt.ebisuzaki > 0
+ inputs1.ebitaus = ebitaus;
+ else
+ inputs1.ebitaus = [];
+ end
+ inputs1.EWSopt = inputs1.EWSopt;
+ inputs1.EWSopt.winsize = sens_opts.winsize;
+ inputs1.EWSopt.bandwidth = sens_opts.bandwidth;
+ inputs1.EWSopt.silent = sens_opts.silent;
+ % res1.indicators = res1.indicators;
+ % res1.colnames = res1.colnames;
+ if nargout > 0
+ res = inputs1;
+ end
+end
+function hfig = plot_sensitiv(results, x, fignr, climit, levels, defaultpos, atitle, EWSopt)
+ hfig = figure;
+ % pcolor(bandwidths, winsizes, x);
+ % shading flat
+ if isfield(EWSopt,'distribution')&&strcmp(EWSopt.distribution, 'uniform')
+ hold on
+ percentiles = [0.05:0.025:0.95];
+ maxcolor = [0.1 0.1 0.1];
+ mincolor = [0.9 0.9 0.9];
+ cs = 0:find(abs(percentiles - 0.5) < 0.01) - 1;
+ cs = cs ./ (max(cs(:)));
+ cs = [cs(1:end - 1), fliplr(cs)];
+ data = quantile(results.indicators(:, :, fignr), percentiles)';
+ ndx = ~isnan(data(:, 1));
+ for i = 1:length(percentiles) - 1
+ c = mincolor + cs(i) .* (maxcolor - mincolor);
+ fill([EWSopt.time(ndx); flipud(EWSopt.time(ndx))], [data(ndx, i); flipud(data(ndx, i + 1))], c);
+ end
+ set(get(gca, 'children'), 'EdgeColor', 'none', 'tag', 'perc_fill')
+ plot(EWSopt.time, quantile(results.indicators(:, :, fignr), 0.5), 'color', maxcolor, 'tag', 'median');
+ xlabel('Time');
+
+ xlims = [min(EWSopt.time), max(EWSopt.time)];
+ xlim(xlims)
+ text(0.1, 0.9, results.colnames{fignr}, 'units', 'normalized', 'fontsize', 11, 'fontname', 'Arial', 'tag', 'toptext')
+ adjustplot('-d');
+ %plot(EWSopt.time,results.indicators(:, :, fignr))
+ else
+ if ~isempty(levels)
+ contourf(results.bandwidths, results.winsizes, x, levels);
+ else
+ contourf(results.bandwidths, results.winsizes, x);
+ end
+ if ~isempty(climit)
+ set(gca, 'clim', climit)
+ end
+ if ~isempty(defaultpos)
+ hold on
+ h1 = plot(defaultpos(1), defaultpos(2), 'k*');
+ set(h1, 'markersize', 6);
+ hold off
+ end
+ adjustplot('-d');
+ xlim([min(results.bandwidths(:)) max(results.bandwidths(:))]);
+ ylim([min(results.winsizes(:)) max(results.winsizes(:))]);
+ xlabel(sprintf('filtering bandwidth (%s)', EWSopt.bandwidth_unit{1}))
+ ylabel(sprintf('sliding window (%s)', EWSopt.winsize_unit));
+ h = colorbar;
+ adjustplot('-d');
+ ylabel(h, atitle);
+ end
+ % title(atitle);
+ set(hfig, 'tag', atitle);
+end
+function pvalues = getpvalues(taus, ebitaus)
+ nebisuzaki = size(ebitaus, 4);
+ pvalues = zeros(size(taus));
+ for i1 = 1:size(taus, 1)
+ for j1 = 1:size(taus, 2)
+ pvalues(i1, j1) = (sum(taus(i1, j1) < ebitaus(i1, j1, 1, :)) + 1) / nebisuzaki;
+ end
+ end
+end
+function EWSopt = mergeoptions(EWSopt, arg)
+ f = fieldnames(EWSopt);
+ if ~isempty(arg) && isstruct(arg{1})
+ %all options can be entered as struct
+ if exist('mergestructs', 'file') == 0
+ initgrind;
+ end
+ %replace common fields in EWSopt
+ EWSopt = mergestructs(EWSopt, arg{1});
+ else
+ for i = 1:2:length(arg)
+ if ~any(strcmp(f, arg{i}))
+ s = sprintf('"%s", ', f{:});
+ error('Unknown option for generic_ews_sens: %s\nValid options are: %s', arg{i}, s(1:end - 1));
+ end
+ if iscell(arg{1})
+ EWSopt.(arg{1}{i}) = arg{1}{i + 1};
+ else
+ EWSopt.(arg{i}) = arg{i + 1};
+ end
+ end
+ end
+end
diff --git a/run_regime.m b/run_regime.m
new file mode 100644
index 0000000..20e6439
--- /dev/null
+++ b/run_regime.m
@@ -0,0 +1,81 @@
+function res=run_regime(tree,regimename)
+ %how to run different regimes
+
+ %regimes, [start year, end year]
+ regimes = {
+ 'BM III', [501 684]; ...
+ 'P I', [705 871]; ...
+ 'P II', [873 1118]; ...
+ 'Early P III', [1157 1215]; ...
+ 'Late P III', [1217 1283]; ...
+ 'P III', [1157 1283]};
+
+ regimendx=find(strcmp(regimes(:,1),regimename));
+ if isempty(regimendx)
+ error('regime "%s" does not exist, check the spelling',regimename);
+ end
+
+ ndx=regimes{regimendx,2};
+ ndx=ndx(1):ndx(2);
+ %ndx= the range of data points for each regime
+ data = [tree.trees(ndx) tree.tot_hab_cell(ndx)];
+
+ if ~any(strcmp(regimename, {'Early P III', 'Late P III'}))
+ %settings for large data sets (all regimes except 'Early PIII' and 'Late PIII'):
+ res = generic_ews(data, 'ebisuzaki', 1000, ...
+ 'detrending' , {'gaussian' 'movmean'}, ...
+ 'bandwidth', {[30] [5 0]}, ...
+ 'winsize', 60, ...
+ 'winsize_unit', 'yr', ...
+ 'bandwidth_unit', {'yr' 'yr'}, ...
+ 'datacolumn', [1 2 1 1 2 2], ...
+ 'arlag', 1, ...
+ 'indicators', {'AR' 'AR' 'crosscorr' 'std' 'std' 'mean'}...
+ );
+
+ res2 = generic_ews_sens(data, 'bandwidth', linspace(5, 40, 15),...
+ 'winsize', linspace(20, 90, 15),...
+ 'winsize_unit', 'yr' ,...
+ 'bandwidth_unit', 'yr',...
+ 'ebisuzaki', 0, ... %put to 1000 for p plot (takes long)
+ 'detrending', 'gaussian', ...
+ 'datacolumn', 1, ...
+ 'arlag', 1, ...
+ 'defaultpos', {[30,60]}, ...
+ 'indicators', {'AR', 'crosscorr', 'std'});
+ figure
+ generic_ews_sens('-p',res2);
+ else
+
+ %settings for small data sets (Early PIII and Late PIII):
+ res = generic_ews(data, 'ebisuzaki', 1000, ...
+ 'detrending' , {'gaussian' ,'movmean'}, ...
+ 'bandwidth', {[15] [5 0]}, ...
+ 'winsize', 20, ...
+ 'winsize_unit', 'yr', ...
+ 'bandwidth_unit', {'yr' 'yr'}, ...
+ 'datacolumn', [1 2 1 1 2 2], ...
+ 'arlag', 1, ...
+ 'indicators', {'AR','AR','crosscorr','std','std','mean'}...
+ );
+
+ res2 = generic_ews_sens(data, 'bandwidth', linspace(5, 40, 15)/2,...
+ 'winsize', linspace(20, 90, 15)/2,...
+ 'winsize_unit', 'yr' ,...
+ 'bandwidth_unit', 'yr',...
+ 'ebisuzaki', 0, ... %put to 1000 for p plot (takes long)
+ 'detrending', 'gaussian', ...
+ 'datacolumn', 1, ...
+ 'arlag', 1, ...
+ 'defaultpos', {[15,20]}, ...
+ 'indicators', {'AR', 'crosscorr', 'std'});
+ figure
+ generic_ews_sens('-p',res2);
+ end
+
+
+
+
+
+end
+
diff --git a/violinplot.m b/violinplot.m
new file mode 100644
index 0000000..cc6ac70
--- /dev/null
+++ b/violinplot.m
@@ -0,0 +1,110 @@
+%violinplot - create a violinplot
+%
+%Usage:
+% violinplot(x,cat) make violinplot of x, where cat are the values of the
+% catagories
+% violinplot(x,cat,'propertery',...)
+% You can use all properties of the function ksdensity (like 'bandwidth') and all
+% properties of patches (like 'facecolor')
+% additionally:
+% 'percentiles' list of percentiles to indicate in the violins for
+% instance: [0.05 0.5 0.95]
+% 'violinwidth' relative width of the violins relative to the smallest gap
+% in the categories, default 2/3
+% 'categorynames' optionnally you can put strings for categories
+%
+
+% (c) Egbert van Nes
+function violinplot(x, cat, varargin)
+ if nargin > 2 && ischar(cat)
+ varargin = [{cat} varargin];
+ cat = [];
+ end
+ if iscell(x)
+ if isempty(cat)
+ cat = 1:numel(x);
+ end
+ x1 = [];
+ cat1 = [];
+ if iscell(cat)
+ for i = 1:numel(x)
+ x1 = [x1; x{i}]; %#ok<AGROW>
+ cat1 = [cat1; cat{i}]; %#ok<AGROW>
+ end
+ else
+ for i = 1:numel(x)
+ x1 = [x1; x{i}]; %#ok<AGROW>
+ cat1 = [cat1; cat(i) + zeros(size(x{i}))]; %#ok<AGROW>
+ end
+ end
+ x = x1;
+ cat = cat1;
+ end
+ ksdensityopts = {'Kernel' , 'Support' , 'Weights', 'Bandwidth', 'Function', 'BoundaryCorrection', 'Censoring' , 'NumPoints'};
+ if nargin == 1 || isempty(cat)
+ cat = repmat(1:size(x, 2), size(x, 1), 1);
+ end
+ if size(x, 1) > 1 && numel(cat) == size(x, 2)
+ cat = repmat(cat, size(x, 1), 1);
+ end
+ percentiles = 0.5;
+ violinwidth = 2 / 3;
+ categorynames = {};
+ used = false(size(varargin));
+ ksdensityndx = false(size(varargin));
+ for i = 2:2:length(varargin)
+ if strcmpi(varargin{i - 1}, 'percentiles')
+ percentiles = varargin{i};
+ used(i - 1:i) = true;
+ elseif strcmpi(varargin{i - 1}, 'violinwidth')
+ violinwidth = varargin{i};
+ used(i - 1:i) = true;
+ elseif strcmpi(varargin{i - 1}, 'categorynames')
+ categorynames = varargin{i};
+ used(i - 1:i) = true;
+ elseif any(strcmpi(varargin{i - 1}, ksdensityopts))
+ ksdensityndx(i - 1:i) = true;
+ end
+ end
+ ksdensityargs = varargin(ksdensityndx);
+ varargin = varargin(~used & ~ksdensityndx);
+ x = x(:);
+ cat = cat(:);
+ ucat = unique(cat);
+ diffcat = min(diff(ucat));
+ if isempty(diffcat)
+ diffcat = 1;
+ end
+ oldhold = ishold;
+ hold on;
+ for i = 1:length(ucat)
+ x1 = x(cat == ucat(i));
+ try
+ [dens, xdens] = ksdensity(x1, ksdensityargs{:});
+ dens = dens / max(dens) * diffcat * violinwidth / 2;
+ xperc = quantile(x1', percentiles);
+ yperc = interp1(xdens, dens, xperc);
+ fill([ucat(i) - dens fliplr(dens + ucat(i))], [xdens fliplr(xdens)], [0.7, 0.7, 1], 'tag', sprintf('violin(%g)', ucat(i)), varargin{:})
+ for j = 1:length(percentiles)
+ if percentiles(j) == 0.5
+ linestyle = 'k-';
+ else
+ linestyle = 'k:';
+ end
+ plot([ucat(i) - yperc(j), ucat(i) + yperc(j)], [xperc(j), xperc(j)], linestyle);
+ end
+ catch err
+ if ~strcmp(err.identifier, 'stats:mvksdensity:NanX')
+ rethrow(err);
+ end
+ end
+ end
+ if ~isempty(categorynames)
+ set(gca, 'xtick', 1:length(categorynames), 'xticklabel', categorynames)
+ end
+ if oldhold
+ hold on
+ else
+ hold off
+ end
+end
--
GitLab