No BSD License  

Highlights from
Tree-Ring MATLAB Toolbox

Tree-Ring MATLAB Toolbox

by

 

Utility functions for analyzing tree-ring data with MATLAB.

slidets
function slidets
% slidets: slide a floating time series along a reference series to find best match
% slidets;
% Last revised 7-24-99
%
% A selected segment of a "floating" time series is slid year-by-year along a reference 
% time series and the correlation coefficient is computed at each iteration. The best 
% match is defined as that with the highest correlation. An optional bootstrap analysis 
% can be used to check how likely this maximum correlation is by chance alone.<P>
%
% Slidets.m was written to help in crossdating remnant-wood ring-width series against 
% very long tree-ring chronologies from a nearby site.<P>
%
% You get to select the segment of the floating series to slide interactively, with 
% guidance from time series plots of the floating series and variance of the floating 
% series.<P>
%
%*** INPUT
%
% Prompts for names of files containing the reference series. 
% Prompt for width of window for sliding segment of floating series. 
% Prompt for ending year of sliding segment.
% Prompt for bootstrap or not
% Prompt for number of bootstrap samles
%
%*** OUTPUT
%
% Screen plots summarizing results.  Key windows are graph of correlation coefficient 
% against year in reference series, and boxplots showing distribution of correlations.  
% Also a boxplot showing distribution of maximum correlations from bootstrapping.
%
%*** REFERENCES -- none
%*** UW FUNCTIONS CALLED 
% corrone
% crn2vec2
% rwchng
%
%*** TOOLBOXES NEEDED
% Signal processing
%
%*** NOTES
%
% To emphasize similarity in year-to-year changes as opposed to trends, both series
% are filtered at the start into scaled first differences.  The series are filtered 
% backwards and forwards with a 9-weight hamming filter to derive smoothed time 
% series.  The first difference of the original time series is divided by the 
% corresponding value of the smoothed series to get the scaled first difference. 
%
% You select the segment of the floating series to be slid along the reference 
% series.  This series is slid one year at a time along the reference series and
% the correlation is computed at each step.  The correlation is plotted as a function 
% of ending year (in the reference series) of the alignment.  A boxplot of the 
% correlations gives some idea of variability of the correlation coefficient for 
% various alignments.  
%
% A plot of variance of the floating series in a sliding window is given to help in 
% choosing the segment of the floating series to test.  The notion is that a period of 
% high variability is a good candidate for getting a definitive match.
%
% Bootstrapped significance level of maximum correlation is available on option.  
% This option takes a while to run.  The bootstrapping can be done once you 
% have found the maximum sliding correlation of the floating segment with the 
% reference series.  A user-specified number of bootstrap samples are drawn from the 
% floating segment (e.g., 1000 samples). Each bootstrap sample is slid along the 
% reference series and the maximum correlation is recorded.  The distribution of 
% maximum bootstrap correlations can then be compared to the maximum correlation 
% found for the original (non-bootstrapped) data. 



% Notation
% variable 1: floating series
% variable 2: reference series
% u full length original series
% v windowed original series
% w windowed scaled first-diff series
% s standard-deviation of rw change, floating series

close all;

% Hard Code
defwidth=50; % default window width
dwin=defwidth; % initialize window width

%--- INPUT THE TIME SERIES

%-- Floating Series -- specify type of file
ksrc=menu('Choose Type of Source File for Floating Series',...
   'Index in ITRDB-formatted .crn file',...
   'Ring width in single .rw file','Two-column time series mtx, year as col 1');
if ksrc==1;
   src = 'crn';
elseif ksrc==2;
   src = 'rw';
elseif ksrc==3; 
   src = 'dat';
else;
   error('ksrc must be 1, 2 or 3');
end
clear ksrc;

%---Floating Series, get the file name
if strcmp(src,'crn');
   [file1,path1]=uigetfile('*.crn','Infile with chronology index');
elseif strcmp(src,'rw');
   [file1,path1]=uigetfile('*.rw','Infile with measured ring width series');
elseif strcmp(src,'dat');
   [file1,path1]=uigetfile('*.dat','Infile of 2-col time series matrix');
end;
pf1=[path1 file1]; % merge path and filename

%--- Floating Series, read the input time series;  
if strcmp(src,'crn');
   [x,s,yr]=crn2vec2(pf1); %x is index, s is sample size, yr is year
elseif strcmp(src,'rw');
   [X,guy,day,pf1]=rwread3(path1,file1); % X has year in col 1, data in col 2
   yr=X(:,1); x=X(:,2);
   clear guy day X;
elseif strcmp(src,'dat');
   eval(['load ' pf1  ' -ascii;']);
   filepref = strtok(file1,'.');
   eval(['x = ' filepref '(:,2);']);
   eval(['yr = ' filepref '(:,1);']);
   clear filepref;
end;

%--- Floating series: store data, name  & year for later use
u1=x; yru1=yr; name1=file1;


%-- Reference Series -- specify type of file
ksrc=menu('Choose Type of Source File for Reference Series',...
   'Index in ITRDB-formatted .crn file',...
   'Ring width in single .rw file','Two-column time series mtx, year as col 1');
if ksrc==1;
   src = 'crn';
elseif ksrc==2;
   src = 'rw';
elseif ksrc==3; 
   src = 'dat';
else;
   error('ksrc must be 1, 2 or 3');
end


%--- Reference Series, get the file name
if strcmp(src,'crn');
   [file2,path2]=uigetfile('*.crn','Infile with chronology index');
elseif strcmp(src,'rw');
   [file2,path2]=uigetfile('*.rw','Infile with measured ring width series');
elseif strcmp(src,'dat');
   [file2,path2]=uigetfile('*.dat','Infile of 2-col time series matrix');
end;
pf2=[path2 file2]; % merge path and filename

%--- Reference Series, read the input time series;  
if strcmp(src,'crn');
   [x,s,yr]=crn2vec2(pf2); %x is index, s is sample size, yr is year
elseif strcmp(src,'rw');
   [X,guy,day,pf2]=rwread3(path2,file2); % X has year in col 1, data in col 2
   yr=X(:,1); x=X(:,2);
   clear guy day X;
elseif strcmp(src,'dat');
   eval(['load ' pf2  ' -ascii;']);
   filepref = strtok(file2,'.');
   eval(['x = ' filepref '(:,2);']);
   eval(['yr = ' filepref '(:,1);']);
   clear filepref;
end;

%--- Reference series: store data, name  & year for later use
u2=x; yru2=yr; name2=file2;

%--- Cleanup
clear x yr src ksrc

%--- COMPUTE SCALED FIRST DIFF SERIES
w1=rwchng(u1);  
w2=rwchng(u2);
yrw1=yru1;  yrw1(1)=[];
yrw2=yru2;  yrw2(1)=[];
nw1=length(w1);
nw2=length(w2);

%--- MAKE WINDOWS WITH TIME SERIES PLOTS OF FULL ORIGINAL SERIES
figure(1); % floating series
plot(yru1,u1); 
title(['Floating Time Series, ' name1]);
xlabel('Year');
grid;
zoom xon;

figure(2); % reference series
plot(yru2,u2); 
title(['Reference Time Series, ' name2]);
xlabel('Year');
grid;
zoom xon;

%--- MAKE WINDOWS WITH FULL LENGTH SCALED FIRST DIFFERENCES
figure(3); % floating series
plot(yrw1,w1); 
title([name1 ': Scaled First-Difference']);
xlabel('Year');
grid;
zoom xon;

figure(4); % reference series
plot(yrw2,w2); 
title([name2 ': Scaled First-Difference']);
xlabel('Year');
grid;
zoom xon;



%--- COMPUTE AND PLOT SLIDING STANDARD DEVIATION OF FLOATING SERIES
% use defwidth-year window

% Compute ending year or each period
iend=nw1:-5:dwin;
iend=(fliplr(iend))'; % to col vector, increasing
nper = length(iend); % number of standard devs
istart = iend-defwidth+1; % start year of period
Wtemp = repmat(NaN,defwidth,nper); % to hold sub-period data
for n = 1:nper;
   i1 = (istart(n):iend(n))';
   Wtemp(:,n)=w1(i1);
end;
s1 = (std(Wtemp))'; % standard deviations for sub-periods
yrs1 = (yrw1(iend));

figure(5);
plot(yrs1,s1);
grid;
title(['Moving Window of Standard Deviation, Scaled First-Difference of ' name1]);
xlabel(['Ending Year of ' int2str(dwin) '-Year Period']);
ylabel('Standard Deviation');
clear Wtemp i1 iend istart nper ;


%--- ANALYSIS LOOP

kwh1 = 1; % while control for level-1 menu
while kwh1==1; 
   
   %--- Set window width
   kmen1 = menu('Choose','Set window width', 'Quit');
   if kmen1 == 1;
      kmen4=menu('Choose',['Use window width of ' int2str(dwin)],'Input new width');
      if kmen4==1; % use current value of window width
         windwid=dwin;
      else;
         prompt={'Enter width'};
         def={int2str(dwin)};
         dlgTitle='Input for width of sliding window (yr)';
         lineNo=1;
         answer=inputdlg(prompt,dlgTitle,lineNo,def);
         windwid=str2num(answer{1});
         dwin=windwid; % re-set default window
      end;
      
      %--- Main analysis 
      
      % Pick test segment
      kmen1a=menu('Choose','Pick end year of test segment','Quit');
      if kmen1a==1; % pick end year of test segment
         kmen1a1=menu('Choose','Graphically--click on end year','In text window');
         if kmen1a1==1; 
            % bring plot of scaled first diffs, floating series, forward
            figure(3);
            xlimits=get(gca,'XLim');
            [t1,w1t1]=ginput(1);
            t1 = round(t1);
            if t1<xlimits(1);
               t1=xlimits(1);
            else;
            end;
            if t1>xlimits(2);
               t1=xlimits(2);
            else;
            end;
            yrend=t1;
            clear t1 xlimits;
         else;
            prompt={'Enter end "year"'};
            def={int2str(max(yru1))}; % default is last year of floating series
            dlgTitle='Input last year of test segment (as numbered in floating series)';
            lineNo=1;
            answer=inputdlg(prompt,dlgTitle,lineNo,def);
            yrend=str2num(answer{1});
         end; % if kmen1a1==1
         
         %--- Check that selected segment valid
         if yrend<yrw1+windwid-1;
            error([int2str(windwid) '-yr period ending in ' int2str(yrend) ' impossible']);
         end;
         yrstart = yrend-windwid+1;
         str1 = sprintf('%5.0f-%5.0f',yrstart,yrend);
         
                 
         kwh2=1; % while control for different test segment
         while kwh2==1;
            % Pull segment of floating series
            L1=yrw1>=yrstart & yrw1<=yrend;
            v1=w1(L1); yrv1=yrw1(L1);
            
            % Make matrix of reference series
            nH=nw2-windwid+1; % number of segments
            Z= repmat(NaN,windwid,nH);
            j1=(0:(windwid-1))';  % column vector
            J1 = repmat(j1,1,nH); 
            j2 = 1:nH; % row vector
            J2 = repmat(j2,windwid,1);
            J3 = J1 + J2; % each col has row indices into Z
            for n = 1:nH;
               j3=J3(:,n);
               Z(:,n)=w2(j3);
            end;
            yrZ = ((yrw2(1)+windwid-1):yrw2(nw2))'; % col vect of ending years of segments
            
            % Compute correlations
            r = corrone(v1,Z);
            
            % Correlations plot
            figure(6);
            plot(yrZ,r);
            grid;
            zoom xon; 
            title(['Correlation of segment ' str1 ' of ' name1]);
            xlabel(['End Year of Segment in ' name2]);
            ylabel('Correlation coefficient');
            
            % Box plot
            figure(7);
            boxplot(r,1,'+',1,1)
            title(['Boxplot of Correlation Coefficients']);
            
            % Top 10 alignments
            figure(8);
            [rsort,isort]=sort(r);
            n10=min([length(isort) 10]);
            rsort=flipud(rsort');
            isort=flipud(isort');
            S1 = ['Correlations with ' name2 ' segments ending in year:'];
            for n = 1:n10;
               str2=sprintf('%5.0f   %6.2f',yrZ(isort(n)), rsort(n));
               S1=strvcat(S1,str2);
            end;
            rwow=rsort(1);
            yrwow=yrZ(isort(1));
            yrbingo=yrZ(isort(1));
            S1=cellstr(S1);
            axes;
            xlims = get(gca,'XLim');
            ylims=get(gca,'YLim');
            xrange=abs(diff(xlims));
            yrange=abs(diff(ylims));
            text(xlims(1)+xrange/10,ylims(2)-yrange/10,S1,...
               'VerticalAlignment','top');
            set(gca,'XTick',[]);
            title([name1 ' Segment ' str1 ' Correlation Summary']);
            
            % Time series plots for period of best match
            figure(9);
            subplot(2,1,1);
            % recall that yrend is the last year in floating segment a
            % recall that yrbingo is the best match year in reference series
            Lu1 = yru1>=yrstart & yru1<=yrend;
            Lw1 = yrw1>=yrstart & yrw1<=yrend;
            yrborn = yrbingo-windwid+1;
            Lu2 = yru2>=yrborn & yru2<=yrbingo;
            Lw2 = yrw2>=yrborn & yrw2<=yrbingo;
            
            % Raw-data plots
            rr=corrcoef([u1(Lu1) u2(Lu2)]);
            rr=rr(1,2);
            strcor1 = sprintf('r = %6.2f',rr);
            plot(yru2(Lu2),zscore(u1(Lu1)),yru2(Lu2),zscore(u2(Lu2)));
            grid;
            title(['Raw Data For Best Match Period: ' strcor1]);
            xlabel('Year');
            ylabel('Z-scores');
            legend(name1,name2);
            zoom xon;
            
            % Scaled First-diff plots
            rr=corrcoef([w1(Lw1) w2(Lw2)]);
            rr=rr(1,2);
            strcor2 = sprintf('r = %6.2f',rr);
            subplot(2,1,2);
            plot(yrw2(Lw2),zscore(w1(Lw1)),yrw2(Lw2),zscore(w2(Lw2)));
            grid;
            title(['Scaled First Differences for Best Match Period: ' strcor2]);
            xlabel('Year');
            ylabel('Z-scores');
            legend(name1,name2);
            zoom xon;
            
            % Bring correlation ts to front
            figure(6);
            
            kmen1a2=menu('Choose','Satisfied','Different test segment?','Quit');
            if kmen1a2==1;
               kmen7=menu('Choose','Bootstrap Analysis','Skip Bootstrapping');
               if kmen7==1; % do bootstrap
                   prompt={'Enter number of samples:'};
                   def={'1000'};
                   dlgTitle='Input desired number of bootstrap samples';
                   lineNo=1;
                   answer=inputdlg(prompt,dlgTitle,lineNo,def);
                   nboot=str2double(answer{1});
                   rboot = repmat(NaN,nboot,1); % to hold bootstrapped r
                   % Generate nboot bootstrap samples of the floating segment
                   % Recall that v1 is the floating segment, with years yrv1
                   [bstat,bsam]=bootstrp(nboot,'mean',v1); % each col of bsam has indices to v1
                   for n = 1:nboot;
                      iboot = bsam(:,n);
                      xboot=v1(iboot);
                      rtemp =max(corrone(xboot,Z)); % highest correl with any segment of ref series
                      rboot(n)=rtemp;
                   end;
                   figure(10);
                   boxplot(rboot);
                   ytemp=get(gca,'YLim');
                   ytemp(2)=max([max(rboot)  min([rwow+.1 1.0])]);
                   set(gca,'YLim',ytemp);
                   text(1,rwow,'o');
                   text(1.1,rwow,'<--- Observed');
                   title(['Maximum r for ' int2str(nboot) ' bootstrapped samples']);
                else; % skip bootstrapping
                end;
             elseif kmen1a2==2; % pick diffent segment
                close(6); close(7); close(8); close(9);
                kwh2=0;
             elseif kmen1a2==3; % Quit 
                kwh2=0;
             end; % kmen1a2==1
          end; % kwh2==1
          
       else;
       end; % if kmen1a==1
       
       
    else;
       kmen1b=menu('Choose','Close windows before quitting','Keep windows open');
       if kmen1b ==1
          close all;
       else;
          % no action- keep windows open
       end;
       kwh1=0;
    end; % if kmen1==1
 end; % while kwh1==1
 

Contact us