clear
close all
figure
set(gcf,'color',[1,1,1])
loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.UEA/.CRU/.TS2p1/.monthly/.mean/.temp/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/T/%28Jan%201981%29%28Dec%202002%29RANGE/T/3/runningAverage/T/%28Jun-Aug%29VALUES/dods')

obs=temp.temp;
loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.EMC/.CFS/.MONTHLY/.surface/.temp/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%28Jan%201981%29%28Dec%202002%29RANGE/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

temp.temp=temp.temp-273.15; 
[nt,m]=size(temp.temp)

h=plot(T/12+1960,obs,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,mean(temp.temp,2),'ko-');
set(h,'MarkerFaceColor','k')
h=plot(T/12+1960,temp.temp,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+prctile(obs,85),'--k');   
hold off
legend('obs','CFS')
title('JJA t2m')

%cfs bias
export_fig -pdf figure1
print -depsc figure1
pause

disp(['correlation =',num2str(corr2(mean(temp.temp,2),obs),2)])

x=[temp.temp,obs];
f=ranks2D(x);
hist(f(:,16),[1:16])
title('rank histogram')
xlabel('rank of the observation in the ensemble')

export_fig -pdf figure1a
print -depsc figure1a
pause
%-------------------------------------------------------------------

cm = mean(temp.temp(:));
co = mean(obs);

h=plot(T/12+1960,obs-co,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,mean(temp.temp,2)-cm,'ko-');
set(h,'MarkerFaceColor','k')
h=plot(T/12+1960,temp.temp-cm,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+prctile(obs-co,85),'--k');   
hold off
legend('obs','CFS')
title('JJA t2m')

%cfs obs anomalies
export_fig -pdf figure2a
print -depsc figure2a
pause

disp(['correlation =',num2str(corr2(mean(temp.temp,2),obs),2)])

x=[temp.temp-cm,obs-co];
f=ranks2D(x);
hist(f(:,16),[1:16])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure2b
print -depsc figure2b
pause


%-------------------------------------------------------------------
h=plot(T/12+1960,obs-co,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,(mean(temp.temp,2)-cm),'ko-');
set(h,'MarkerFaceColor','k')
wrk = mean(temp.temp,2);
wrk2 = (temp.temp -wrk(:,ones(1,15)))/2+wrk(:,ones(1,15));
h=plot(T/12+1960,wrk2-cm,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+prctile(obs-co,85),'--k');   
hold off
legend('obs','CFS/2')
title('JJA t2m')

%cfs obs anomalies spread too small by a factor of 2
export_fig -pdf figure2d
print -depsc figure2d
pause


disp(['correlation =',num2str(corr2(mean(temp.temp,2),obs),2)])

%cfs obs anomalies spread too small by a factor of 2
x=[wrk2-cm,obs-co];
f=ranks2D(x);
hist(f(:,16),[1:16])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure2e
print -depsc figure2e
pause


%-------------------------------------------------------------------
obs_rank = ranks(obs);
obs_prctile = obs_rank/length(obs_rank)*100;

model_rank = ranks(temp.temp(:));
model_prctile = model_rank/length(model_rank)*100;
model_prctile = reshape(model_prctile,nt,m);

mean_rank = ranks([mean(temp.temp,2);temp.temp(:)]); 
mean_prctile=mean_rank(1:nt)/(nt+m*nt)*100;

h=plot(T/12+1960,obs_prctile,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,mean_prctile,'ko-');
set(h,'MarkerFaceColor','k')
h=plot(T/12+1960,model_prctile,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+85,'--k');   
hold off
legend('obs','CFS')
title('JJA t2m percentiles')

%cfs bias fixed by percentiles
export_fig -pdf figure2
print -depsc figure2
pause

x=[model_prctile,obs_prctile'];
f=ranks2D(x);
hist(f(:,16),[1:16])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure2c
print -depsc figure2c
pause
%-------------------------------------------------------------------

loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.UEA/.CRU/.TS2p1/.monthly/.mean/.temp/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/T/%28Jan%201960%29%28Dec%202001%29RANGE/T/3/runningAverage/T/%28Jun-Aug%29VALUES/dods')

obs=temp.temp;

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.ECMWF/.t2m/S/%28Jan%201960%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

mod1=t2m.t2m-273.15; 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.UKMO/.t2m/S/%28Jan%201960%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

mod2=t2m.t2m-273.15; 

h=plot(T/12+1960,obs,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obs,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

% two models with different bias
export_fig -pdf figure3
print -depsc figure3
pause

x=[mod1,mod2,obs];
disp(['correlation=',num2str(corr2(mean(x(:,1:(end-1)),2),x(:,end)),2)])
disp(['model 1 correlation=',num2str(corr2(mean(x(:,1:9),2),x(:,end)),2)])
disp(['model 2 correlation=',num2str(corr2(mean(x(:,10:(end-1)),2),x(:,end)),2)])
f=ranks2D(x);
nn=size(x,2);
hist(f(:,nn),[1:nn])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure3b
print -depsc figure3b

%----------------------------------------------------------------------------
mod1a = mod1-mean([mod1(:);mod2(:)]);
mod2a = mod2-mean([mod2(:):mod2(:)]);
obsa = obs-mean(obs(:));

h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

% two models with different bias single bias removed
export_fig -pdf figure3a
print -depsc figure3a
pause

x=[mod1a,mod2a,obsa];
disp(['correlation=',num2str(corr2(mean(x(:,1:(end-1)),2),x(:,end)),2)])
f=ranks2D(x);
nn=size(x,2);
hist(f(:,nn),[1:nn])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure3c
print -depsc figure3c


%----------------------------------------------------------------------------

mod1a = mod1-mean(mod1(:));
mod2a = mod2-mean(mod2(:));
obsa = obs-mean(obs(:));
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

% two models with different bias individual bias removed
export_fig -pdf figure4
print -depsc figure4
pause

x=[mod1a,mod2a,obsa];
disp(['correlation=',num2str(corr2(mean(x(:,1:(end-1)),2),x(:,end)),2)])
f=ranks2D(x);
nn=size(x,2);
hist(f(:,nn),[1:nn])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure4a
print -depsc figure4a

%----------------------------------------------------------------------------
mod2a = 2*mod2a;
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','2 x model2','location','northwest')
title('JJA t2m')

export_fig -pdf figure5
print -depsc figure5
pause

x=[mod1a,mod2a,obsa];
disp(['correlation=',num2str(corr2(mean(x(:,1:(end-1)),2),x(:,end)),2)])
f=ranks2D(x);
nn=size(x,2);
hist(f(:,nn),[1:nn])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure5a
print -depsc figure5a


%----------------------------------------------------------------------------
mod1a = mod1a/std(mod1a(:));
mod2a = mod2a/std(mod2a(:));
obsa = obsa/std(obsa);
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

export_fig -pdf figure6
print -depsc figure6
pause

x=[mod1a,mod2a,obsa];
disp(['correlation=',num2str(corr2(mean(x(:,1:(end-1)),2),x(:,end)),2)])
f=ranks2D(x);
nn=size(x,2);
hist(f(:,nn),[1:nn])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure6a
print -depsc figure6a
 
%------------------------------------------------------------------------
mod1a = mod1-mean(mod1(:));
mod2a = mod2-mean(mod2(:));
obsa = obs-mean(obs(:));

x=[mod1a,mod2a,obsa];

p=mean(x(:,1:(end-1))>prctile(obsa,85),2);
po=x(:,end)>prctile(obsa,85);

bar(T/12+1960,p)
hold on
plot(T/12+1960,po,'rx')
plot(T(13)/12+1960,po(13),'r*')
hold off
title('Probability of exceeding of 85-th percentile')
legend('forecast prob','observation')

disp(['BSS=',num2str(1- (sum(p-po).^2)./sum((0.15-po).^2),2)])

export_fig -pdf figure7
print -depsc figure7
 
m = mean(x(:,1:(end-1)),2);
s = std(x(:,1:(end-1)),[],2);
pg = 1 - normcdf(prctile(obsa,85),m,s);

bar(T/12+1960,pg)
hold on
plot(T/12+1960,po,'rx')
plot(T(13)/12+1960,po(13),'r*')
hold off
title('Probability of exceeding of 85-th percentile')
legend('forecast prob (Gaussian)','observation')

disp(['BSS=',num2str(1- (sum(pg-po).^2)./sum((0.15-po).^2),2)])

export_fig -pdf figure7a
print -depsc figure7a
 
%constant ensemble
x=[mod1a,mod2a,obsa];
wrk = x(:,1:(end-1)) - m(:,ones(1,18));
wrk = wrk(:)';
ec = m(:,ones(1,length(wrk))) + wrk(ones(length(m),1),:);

h=plot(T/12+1960,ec,'+');   
set(h,'color',[1,1,1]*0.5)

hold on
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')

h=plot(T/12+1960,mean(ec,2),'ko-');
set(h,'MarkerFaceColor','k')
plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','ens')
title('JJA t2m')

export_fig -pdf figure8
print -depsc figure8
pause 

x=[ec,obsa];
disp(['correlation=',num2str(corr2(mean(x(:,1:(end-1)),2),x(:,end)),2)])
f=ranks2D(x)/42;
nn=size(x,2);
hist(f(:,nn),[1:18])
title('rank histogram')
xlabel('rank of the observation in the ensemble')
export_fig -pdf figure8a
print -depsc figure8a
 
p=mean(ec(:,1:(end-1))>prctile(obsa,85),2);

bar(T/12+1960,p)
hold on
plot(T/12+1960,po,'rx')
plot(T(13)/12+1960,po(13),'r*')
hold off
title('Probability of exceeding of 85-th percentile')
legend('forecast prob const. ens','observation')

disp(['BSS=',num2str(1- (sum(p-po).^2)./sum((0.15-po).^2),2)])

export_fig -pdf figure8b
print -depsc figure8b

x=[mod1a,mod2a,obsa];
p=mean(x(:,1:(end-1))>prctile(obsa,85),2);

b=glmfit(m,p,'binomial','probit');
pglm = glmval(b,m,'probit');

bar(T/12+1960,pglm)
hold on
plot(T/12+1960,po,'rx')
plot(T(13)/12+1960,po(13),'r*')
hold off
title('Probability of exceeding of 85-th percentile')
legend('forecast prob glm','observation')

disp(['BSS=',num2str(1- (sum(pglm-po).^2)./sum((0.15-po).^2),2)])

export_fig -pdf figure9
print -depsc figure9
 
mm=-1:0.01:1;
pglm = glmval(b,mm,'probit');
h=plot(mm,pglm);
set(h,'linewidth',1.5)
xlabel('ensemble mean')
ylabel('P(exceeding 85th %ile')
hold on
h=plot(m,p,'ko');
set(h,'markersize',7,'markerfacecolor','k')
h=plot(m,po,'ro');
set(h,'markersize',4,'markerfacecolor','r')

hold off
legend('GLM','ens prob','obs','location','northwest')

export_fig -pdf figure9a
print -depsc figure9a


return


mc=10000;
m=zeros(1,mc);
wrk = temp.temp(~isnan(temp.temp));
for imc = 1:mc
  ind=randperm(length(wrk));
  m(imc)=mean(wrk(ind(1:15)));
end
ml = prctile(m,2.5);
mu = prctile(m,97.5);


plot(S/12+1960,temp.temp,'+');   
hold on
h=plot(S/12+1960,mean(temp.temp,2),'o');
set(h,'MarkerFaceColor','k')
plot(S/12+1960,S*0+ml,'--')
plot(S/12+1960,S*0+mu,'--')
h=plot(T/12+1960,obs,'ro-');
set(h,'MarkerFaceColor','r')
h=plot(S/12+1960,mean(temp.temp,2)-mean(wrk)+mean(obs),'go-');
set(h,'MarkerFaceColor','g')
hold off

corr2(obs,mean(temp.temp(1:22,:),2))

 