crospull.m

function IX=crospull(nyears,nlag,plag)
% crosspull: make logical pointer to estimation and crossvalidation subsets for reconstruction model
% IX=crospull(nyears,nlag,plag);
% Last revised 7-13-99
%
% For a data matrix for climatic reconstruction with a model that possibly includes 
% lags on predictors, makes a logical pointer matrix to subsets of rows (years) to 
% be used in estimation and crossvalidation.  Uses "leave-n-out" strategy, where n 
% is two times the maximum of number of negative and positive lags for the reconstruction
% model.  
%
%*** INPUT ARGUMENTS
%
%  nyears (1x1)i  number of years in the calibration/verification period
%  nlag (1x1)i    number of negative lags on the predictors
%  plag (1x1)i    number of positive lags on the predictors
%
%
%*** OUTPUT ARGUMENTS
%
%  IX (m1 x n1)L   logical array, elements of a col refer to rows of data
%		matrix in the calling program. First column of IX has 1's for rows to 
%     be used 
%
%*** REFERENCES
%
% Crossvalidation is discussed by Michaelsen, J., 1987, Cross-validation in statistical 
% climate forecast models, J. of Climate and Applied Meterorology 26, 1589-1600.
%
% The "leave-n-out" strategy is used in:
% Meko, D.M., 1997, Dendroclimatic reconstruction with time varying subsets of tree 
% indices: Journal of Climate, v. 10, p. 687-696.
%
%*** UW FUNCTIONS CALLED -- none
%*** TOOLBOXES NEEDE -- none
%
%*** NOTES
%
% Makes an intermediate array of indexes to be used in cross
% validation in regression models.  Assume you are doing a regression
% analysis with crossvalidation by a "leave n out" strategy.  There
% are n years of data, and you will repetitively estimate a regression
% equation, leaving out enough years so that the predictand for the
% calibration period is independent of tree-ring data that might
% be related to the key year of the left-out sequence.
%
% Given the number of years of data you have, and the number of positive
% and negative lags on the predictor(s) in the model, this function builds
% a 1/0 logical pointer to your data telling which years to use in the 
% current model, and a subscript array telling which year is to be 
% estimated using the current model.
%
% Calling program would use column one of IX for first model:  rows with 1 would
% be used to fit the model, rows would zero would not be used for the fit, and
% the crossvalidation estimate is for year 1 of the nyears.  Next the calling program  
% would use column 2 of IX:  rows with 1 would be used for the fit, rows 
% with zero would not, and the crossvalidation estimate is for year 2. And so on 
% through all nyears of the data.
%
% It is necessary to insure that the data 
% for the central year of the "left-out" segment is independent of the 
% data used to estimate the parameters.  crospull.m handles this problem
% conservatively by assuming that 2 times the larger of plag and nlag  
% should be the number of years buffer on either side of the central year.
% Thus a model with 1 positive lag and 2 negative lags would look like a 
% model with 2 positive lags and two negative lags.
%
%
%  m1 equals nyears
%  n1 is the number of independent predictions that can be done in cross
%	validation.  n1=m1


maxlag=max(nlag,plag);
nzero = 1 + 4*(maxlag);  %  This number of years will be left out
%	of each regression (except a few) in the crossvalidation procedure.
%  The central
%	year of the left-out set will be the year predicted and checked
%	in crossvalidation.  The exceptions are the first and last 
%   2*maxlag data points, for which fewer than nzero need to be zeroed
%   out.

nmods=nyears;  % Number of regression models to estimate
n1= nmods;

IX = ones(nyears,nmods);  % Initialize logical array as all ones


%*****  Fill an intermediate subscript array

L1=(1:nzero)';  % cv like [1 2 3 4 5]'
L2=L1(:,ones(nmods-4*maxlag,1));  % dupe that column multiple times
K1= 0:(nmods-4*maxlag-1);  %  rv like [0 1 2 3 ...  45]
K2=  K1(ones(nzero,1),:);  % dupe that row multiple times
LK=L2+K2;  % For each regression model, a col of LK tells which rows are
%	to be left out.
%*******  Loop for "inside" years 

first=2*maxlag+1;
last=nyears-(2*maxlag);
for j=1:(nyears-4*maxlag);
	ithis=LK(:,j);
	IX(ithis,j+2*maxlag)=(zeros(nzero,1));
end



%***********   HANDLE "OUTER" YEARS

ncols=2*maxlag;  % this many columns are "outer" on each end
if ncols~=0;  % no outer-year problems if no lags
   xl=ones(nyears,ncols);
   for i=1:ncols;
      for j=1:(2*maxlag+i);
         xl(j,i)=0;
      end
   end
   xr=fliplr(xl);
   xr=flipud(xr);
   IX(:,1:(2*maxlag))=xl;
   IX(:,(nyears-2*maxlag+1):nyears)=xr;
end;

IX=logical(IX);