loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.CPC/.CAMS/.station/.temperature/IWMO/%2816110.%29VALUE/.temp/T/3/runningAverage/T/%281980%29%282009%29RANGE/T/%28Jun-Aug%29VALUES/dods')

t = temp.temp;

[m,s] =normfit(t);

tt = [(m-4*s):0.01:(m+4*s)];

h=plot(tt,normpdf(tt,m,s));
set(h,'linewidth',5);
set(h(1),'color',[0 0 0]+0.3)

hold on
h=plot(t,0*t,'r+');
set(h,'markersize',10)
plot([1,1]*norminv(1/3)*s+m,[0,0.4])
plot([1,1]*norminv(2/3)*s+m,[0,0.4])
hold off
title('JJA temperature in Trieste 1980-2009')
grid
export_fig('-pdf',['trieste_temp'])

% cdf
h=plot(tt,normcdf(tt,m,s));
set(h,'linewidth',5);
set(h(1),'color',[0 0 0]+0.3)
grid
hold on
h=plot(t,0*t,'r+');
set(h,'markersize',10)
plot([[1,1]*norminv(1/3)*s+m,tt(1)-1],[0,1/3,1/3])
plot([[1,1]*norminv(2/3)*s+m,tt(1)-1],[0,2/3,2/3])
hold off
set(gca,'xlim',[tt(1),tt(end)])
title('JJA temperature in Trieste CDF (1980-2009)')
grid
export_fig('-pdf',['trieste_temp_cdf'])

%modest skill forecast
h=plot(tt,normpdf(tt,m,s));
set(h,'linewidth',5);
set(h(1),'color',[0 0 0]+0.3)

hold on
r=0.4;
sig = r*s;
noise = sqrt(1-r^2)*s;
h=plot(tt,normpdf(tt,m+sig,noise));
set(h,'linewidth',5);
set(h(1),'color',[1 0 0])
h=plot(t,0*t,'r+');
set(h,'markersize',10)
plot([1,1]*norminv(1/3)*s+m,[0,0.4])
plot([1,1]*norminv(2/3)*s+m,[0,0.4])
hold off
title('JJA temperature in Trieste 1980-2009')
grid
export_fig('-pdf',['trieste_temp_forecast1'])

100*[normcdf(norminv(1/3)*s+m,m+sig,noise),...
1-normcdf(norminv(2/3)*s+m,m+sig,noise),...
normcdf(norminv(2/3)*s+m,m+sig,noise)-normcdf(norminv(1/3)*s+m,m+sig,noise)]


%high skill forecast
h=plot(tt,normpdf(tt,m,s));
set(h,'linewidth',5);
set(h(1),'color',[0 0 0]+0.3)

hold on
r=0.8;
sig = r*s;
noise = sqrt(1-r^2)*s;
h=plot(tt,normpdf(tt,m+sig,noise));
set(h,'linewidth',5);
set(h(1),'color',[1 0 0])

h=plot(t,0*t,'r+');
set(h,'markersize',10)
plot([1,1]*norminv(1/3)*s+m,[0,0.7])
plot([1,1]*norminv(2/3)*s+m,[0,0.7])
hold off
title('JJA temperature in Trieste 1980-2009')
grid
export_fig('-pdf',['trieste_temp_forecast2'])

100*[normcdf(norminv(1/3)*s+m,m+sig,noise),...
1-normcdf(norminv(2/3)*s+m,m+sig,noise),...
normcdf(norminv(2/3)*s+m,m+sig,noise)-normcdf(norminv(1/3)*s+m,m+sig,noise)]
