Documentation of regress

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

clear

Cross-Reference Information

Listing of script regress

cd /home/disk/hayes2/dvimont/csm/300yr/srf
load sst1-150.mat
load oro1-150.mat
load lims1-150.mat
whos

lims = [0 360 -90 90];
%sst = sst(360:1800,:,:);
%oro = oro(360:1800,:,:);
sst = subset(sst, nxa, nya, lims);
oro = subset(oro, nxa, nya, lims);
nxkeep = nxa(find(nxa >= lims(1) & nxa <= lims(2)));
nykeep = nya(find(nya >= lims(3) & nya <= lims(4)));
[ntmstp, nypts, nxpts] = size(sst);
whos

sst = reshape(sst, ntmstp, nypts*nxpts);
oro = reshape(oro, ntmstp, nypts*nxpts);
oro = mean(oro);
dokeep = find(oro == 0);
sst = sst(:,dokeep);

%
%  Alter sst for ice points.  Notice, the annual cycle feature was
%    not included for figures 1 and 2 in ~dvimont/pacific
%

%noelim = [];
%for i = 1:length(dokeep);
%  noelimtmp = [];
%  for j = 1:12
%    icept = (find(sst(j:12:ntmstp,i) <= 273.14)-1).*12 + j;
%    noice = (find(sst(j:12:ntmstp,i) >= 273.14)-1).*12 + j;
%    if (~isempty(icept) & ~isempty(noice));
%      sst(icept,i) = mean(sst(noice,i)) * ones(length(icept),1);
%    end
%    if ~isempty(noice);
%      noelimtmp = [noelimtmp j];
%    else
%      sst(icept,i) = mean(sst(icept,i)) * ones(length(icept),1);
%    end
%  end
%  if ~isempty(noelimtmp);
%    noelim = [noelim i];
%  end
%end
%sst = sst(:,noelim);
%dokeep = dokeep(noelim);

noelim = [];
for i = 1:length(dokeep);
  icept = find(sst(1:ntmstp,i) <= 273.14);
  if isempty(icept);
    noelim = [noelim i];
  end
end
sst = sst(:,noelim);
dokeep = dokeep(noelim);
whos

%
%  Remove linear trend from each point's time series
%

x = 1:ntmstp;
xprime = x - mean(x);
a1 = xprime * (sst - ones(ntmstp,1)*mean(sst)) / (xprime * xprime');
a0 = mean(sst) - a1*mean(x);
yhat = ones(ntmstp,1)*a0 + x'*a1;

sst = sst - yhat;
%yhat = yhat - ones(ntmstp,1)*mean(yhat);
%for i=1:length(dokeep)
%  yhat(:,i) = yhat(:,i) / sqrt((yhat(:,i))'*yhat(:,i));
%end

clear yhat;

%  Remove climatology

[sst, clim] = annave(sst);

%trendpat = zeros(1,length(dokeep));
%for i=1:length(dokeep);
%  trendpat(1,i) = (yhat(:,i))' * sst(:,i);
%end

%temp = zeros(1, nxpts*nypts);
%temp(1,dokeep) = trendpat;
%gcont(squeeze(reshape(temp,1,nypts,nxpts)),lims,10)

sstemp = zeros(ntmstp, nxpts*nypts);
sstemp(:,dokeep) = sst;
sst = sstemp;
clear sstemp

sst = cosweight(sst, nya);
sst = sst(:, dokeep);

%  Running Mean

sst3yr = myrunning_ave(sst, 36);

%  This might make things faster

sst3yr = sst3yr(1:3:ntmstp,:);

%  Take EOF's.  (Might alter the dimensions if need be...)

c = sst3yr * sst3yr' / (length(dokeep));
[lam, eof, per] = eof(c);

pc1 = eof(:,1:10);
eof1 = pc1'*sst3yr;
eof1 = eof1' * diag(1 ./ sqrt(diag(eof1*eof1'/length(dokeep))));
eof1 = eof1';

rescale = std(pc1);
pc1 = pc1 * diag(1./rescale);
eof1 = diag(rescale) * eof1;

temp = zeros(10, nypts*nxpts);
temp(1:10,dokeep) = eof1(1:10,:);

%
%  Determine significance of eigenvalues
%

for i = 1:length(dokeep);
  a(i) = 0.5*(corr(sst3yr(:,i),sst3yr(:,i),1) + sqrt(corr(sst3yr(:,i),sst3yr(:,i),2)));
end
a = mean(a);
dof = -1 * 600 * log(a) / 2;

dlam = lam * sqrt(2/dof);