clear
loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/.monthly/.sst/T/3/runningAverage/T/%28Dec-Feb%29VALUES/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/.climatology/.mask/1/masklt/add/Y/-20/20/RANGEEDGES/dods')

sst = asum.asum;

y = squeeze(sst(Y==-11.5,X==195.5,:));

[ny,nx,nt] = size(sst);
wrk = reshape(sst,ny*nx,nt);
wm = mean(wrk,2);
wrk = wrk - wm(:,ones(1,nt));
mask = ~isnan(wrk(:,1));
[uu,ss,vv] = svd(wrk(mask,:),'econ');

[corr2(vv(:,1),y),corr2(vv(:,2),y),corr2(vv(:,3),y),corr2(vv(:,4),y)]

for im=1:100
%generate 30 years of data
nt=30

fake = uu*ss*randn(29,nt) + wm(mask,ones(1,nt));
wrk=zeros(ny*nx,nt)+NaN;
wrk(mask,:)=fake;
sst = reshape(wrk,ny,nx,nt);

y = squeeze(sst(Y==-11.5,X==195.5,:));


[ny,nx,nt] = size(sst);
wrk = reshape(sst,ny*nx,nt);
wm = mean(wrk,2);
wrk = wrk - wm(:,ones(1,nt));
mask = ~isnan(wrk(:,1));
[u,s,v] = svd(wrk(mask,:),'econ');
[corr2(v(:,1),y),corr2(v(:,2),y),corr2(v(:,3),y),corr2(v(:,4),y)]

XX = [v(:,1),ones(nt,1)];
V = XX*inv(XX'*XX)*XX';%HAT matrix
res = (eye(nt)-V)*y;
res_cv = res./(1-diag(V));%cv residual w/o cv!!!!!!!!!!!!!! see book
b = regress(y,XX);

cvRes(im) = mean(res_cv.^2);


%generate 1000 years of data
nt=1000;

fake = uu*ss*randn(29,nt) + wm(mask,ones(1,nt));
wrk=zeros(ny*nx,nt)+NaN;
wrk(mask,:)=fake;
sst = reshape(wrk,ny,nx,nt);

y = squeeze(sst(Y==-11.5,X==195.5,:));

yp = [(u(:,1)'*(fake-wm(mask,ones(1,nt)))/s(1,1))',ones(nt,1)]*b;

Res(im) = mean((y-yp).^2);

mean(cvRes<Res)
end