clear all
%
fs1 = 15;

fs2 = 12;

nx = 500;
x = linspace(0,1,nx);
y1 = linspace(0,1,nx);
y = fliplr(y1);
%
l2=log10(2.0)

load nw.tmp;

ar = nw(:, 6:6);
ca = nw(:, 7:7);
nu = nw(:, 8:8);

nnw=sum(nu)

arms=sum(ar)
cams=sum(ca)
ar1 = ar/arms;
ca1 = ca/cams;
arm = cumsum(ar)/arms;
cam = cumsum(ca)/cams;
%
de2 = ca1./ar1;
lde2nw = log10(de2)/l2;
smlde2nw = sum(lde2nw.*ca1)
%
smlde2scnw = sqrt(sum((lde2nw-smlde2nw).^2.*ca1))
smlde2scnw3 = sum((lde2nw-smlde2nw).^3.*ca1)
smlde2scnw4 = sum((lde2nw-smlde2nw).^4.*ca1)
snw = smlde2scnw/sqrt(nnw)
skewnw=smlde2scnw3/smlde2scnw^3
kurtnw=smlde2scnw4/smlde2scnw^4-3
%
nu0=0.5
inf=2.3647
tinf=2^inf
%
% d=2;
% d=3;
% d=4;
% d=5;
% d=10;
% d=20;
% d=50;
% options=optimset('Display','iter');
%
% [nu,fval] = fsolve(@(nu) twop(nu,d),nu0,options)
%
d=tinf;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu1=[1, 1-nu, 0];
tau1=[0, nu/d, 1];
nus=nu1;
taus=tau1;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(1)=d;
% sig(1)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(1)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(1)=0
kurt(1)=0

d=tinf*1.10;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu2=[1, 1-nu, 0];
tau2=[0, nu/d, 1];
nus=nu2;
taus=tau2;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(2)=d;
% sig(2)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(2)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(2)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(2)^3
kurt(2)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(2)^4-3

d=tinf*1.25;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu3=[1, 1-nu, 0];
tau3=[0, nu/d, 1];
nus=nu3;
taus=tau3;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(3)=d;
% sig(3)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(3)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(3)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(3)^3
kurt(3)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(3)^4-3

d=tinf*1.5
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu4=[1, 1-nu, 0];
tau4=[0, nu/d, 1];
nus=nu4;
taus=tau4;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(4)=d;
% sig(4)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(4)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(4)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(4)^3
kurt(4)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(4)^4-3

d=tinf*2
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu5=[1, 1-nu, 0];
tau5=[0, nu/d, 1];
nus=nu5;
taus=tau5;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(5)=d;
% sig(5)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(5)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(5)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(5)^3
kurt(5)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(5)^4-3

d=tinf*3;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu6=[1, 1-nu, 0];
tau6=[0, nu/d, 1];
nus=nu6;
taus=tau6;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(6)=d;
% sig(6)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(6)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(6)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(6)^3
kurt(6)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(6)^4-3

d=tinf*5;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu7=[1, 1-nu, 0];
tau7=[0, nu/d, 1];
nus=nu7;
taus=tau7;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(7)=d;
% sig(7)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(7)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(7)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(7)^3
kurt(7)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(7)^4-3

d=tinf*10;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu8=[1, 1-nu, 0];
tau8=[0, nu/d, 1];
nus=nu8;
taus=tau8;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(8)=d;
% sig(8)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(8)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(8)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(8)^3
kurt(8)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(8)^4-3

d=tinf*50;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu9=[1, 1-nu, 0];
tau9=[0, nu/d, 1];
nus=nu9;
taus=tau9;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(9)=d;
% sig(9)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(9)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(9)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(9)^3
kurt(9)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(9)^4-3

d=tinf*100;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu10=[1, 1-nu, 0];
tau10=[0, nu/d, 1];
nus=nu10;
taus=tau10;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(10)=d;
% sig(10)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(10)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(10)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(10)^3
kurt(10)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(10)^4-3

d=tinf*250;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu11=[1, 1-nu, 0];
tau11=[0, nu/d, 1];
nus=nu11;
taus=tau11;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(11)=d;
% sig(11)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(11)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(11)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(11)^3
kurt(11)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(11)^4-3

d=tinf*1000;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu12=[1, 1-nu, 0];
tau12=[0, nu/d, 1];
nus=nu12;
taus=tau12;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(12)=d;
% sig(12)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(12)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(12)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(12)^3
kurt(12)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(12)^4-3

d=tinf*10000;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu13=[1, 1-nu, 0];
tau13=[0, nu/d, 1];
nus=nu13;
taus=tau13;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(13)=d;
% sig(13)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(13)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(13)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(13)^3
kurt(13)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(13)^4-3

d=tinf*100000;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu14=[1, 1-nu, 0];
tau14=[0, nu/d, 1];
nus=nu14;
taus=tau14;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(14)=d;
% sig(14)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(14)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(14)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(14)^3
kurt(14)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(14)^4-3

d=tinf*1000000;
[nu,fval] = fsolve(@(nu) twop1(nu,d,inf),nu0);
nu15=[1, 1-nu, 0];
tau15=[0, nu/d, 1];
nus=nu15;
taus=tau15;
y1=log10(d)/l2
y2=log10(nus(2)/(1-taus(2)))/l2
infs=y1*(1-nus(2))+y2*nus(2)
dr(15)=d;
% sig(15)=sqrt(y1^2*(1-nus(2))+y2^2*nus(2)-infs^2)
sig(15)=sqrt((y1-infs)^2*(1-nus(2))+(y2-infs)^2*nus(2))
skew(15)=((y1-infs)^3*(1-nus(2))+(y2-infs)^3*nus(2))/sig(15)^3
kurt(15)=((y1-infs)^4*(1-nus(2))+(y2-infs)^4*nus(2))/sig(15)^4-3

%
figure
% plot(tau2, nu2)
plot(x,y,'k-',tau1,nu1,'m:',tau2,nu2,'b:',tau3,nu3,'g:',tau4,nu4,'c:',tau5,nu5,'m--',....
tau6,nu6,'b--',tau7,nu7,'g--',tau8,nu8,'c--',tau9,nu9,'m-',tau10,nu10,'b-',tau11,nu11,'g-',...
tau12,nu12,'c-', tau13,nu13,'m:',tau14,nu14,'b:',tau15,nu15,'g:', arm,1-cam,'r-', 'LineWidth', 2.00)
%
grid
xtxt = 0.1
ytxt = 0.95
st=0.1
%
% text (xtxt, ytxt, ['  NW info (r - bits) = ' num2str(smlde2nw,6),...
% ',  \sigma = ' num2str(smlde2scnw,4)], 'FontSize', fs2)
% text (xtxt, ytxt-st, [' N = ' num2str(nnw,4),...
% ', \sigma N^{-0.5} = ' num2str(snw,4)], 'FontSize', fs2)               %
% 
axis square

title('Fig. 8', 'FontSize', fs1)
%title('Error diagram,  NW Pacific, info content -- 2.36 bits', 'FontSize', fs1)
%
xlabel('Fraction alarm area, \tau', 'FontSize', fs1)
% ylabel('\nu', 'FontSize', fs1)
ylabel('Fraction failures to predict, \nu', 'FontSize', fs1)
figure
% plot(tau2, nu2)
plot(x,y,'k-',tau1,nu1,'m:',tau2,nu2,'b:',tau3,nu3,'g:',tau4,nu4,'c:',tau5,nu5,'m--',....
tau6,nu6,'b--',tau7,nu7,'g--',tau8,nu8,'c--',tau9,nu9,'m-',tau10,nu10,'b-',tau11,nu11,'g-',...
tau12,nu12,'c-', tau13,nu13,'m:',tau14,nu14,'b:',tau15,nu15,'g:', arm,1-cam,'r-', 'LineWidth', 2.00)

%
grid
xtxt = 0.01
ytxt = 0.975
st=0.05
%
text (xtxt, ytxt, ['  NW info (r - bits) = ' num2str(smlde2nw,6),...
',  \sigma = ' num2str(smlde2scnw,4)], 'FontSize', fs2)
text (xtxt, ytxt-st, [', N = ' num2str(nnw,4),...
', \sigma N^{-0.5} = ' num2str(snw,4)], 'FontSize', fs2)               %

axis square
axis([0 0.1 0.4 1])
%
title('Error diagram,  NW Pacific, info content -- 2.36 bits', 'FontSize', fs1)

xlabel('Fraction alarm area, \tau', 'FontSize', fs1)
% ylabel('\nu', 'FontSize', fs1)
ylabel('Fraction failures to predict, \nu', 'FontSize', fs1)


figure
% plot(tau2, nu2)
plot(x,y,'k-',tau1,nu1,'m:',tau2,nu2,'b:',tau3,nu3,'g:',tau4,nu4,'c:',tau5,nu5,'m--',....
tau6,nu6,'b--',tau7,nu7,'g--',tau8,nu8,'c--',tau9,nu9,'m-',tau10,nu10,'b-',tau11,nu11,'g-',...
tau12,nu12,'c-', tau13,nu13,'m:',tau14,nu14,'b:',tau15,nu15,'g:', arm,1-cam,'r-', 'LineWidth', 2.00)

%
grid
xtxt = 0.001
ytxt = 0.985
st=0.02
%
text (xtxt, ytxt, ['  NW info (r - bits) = ' num2str(smlde2nw,6),...
',  \sigma = ' num2str(smlde2scnw,4)], 'FontSize', fs2)
text (xtxt, ytxt-st, [', N = ' num2str(nnw,4),...
', \sigma N^{-0.5} = ' num2str(snw,4)], 'FontSize', fs2)               %

axis square
axis([0 0.01 0.75 1])
%
title('Error diagram,  NW Pacific, info content -- 2.36 bits', 'FontSize', fs1)

%
xlabel('Fraction alarm area, \tau', 'FontSize', fs1)
% ylabel('\nu', 'FontSize', fs1)
ylabel('Fraction failures to predict, \nu', 'FontSize', fs1)


figure
xl=10000000
x3=[1,xl];
y2a=[smlde2scnw, smlde2scnw];
y3=[skewnw, skewnw];
y4=[kurtnw, kurtnw];

semilogx (dr,sig,'r-o', dr,skew,'b-o', dr,kurt,'g-o', x3,y2a,'k-',...
x3,y3,'b-', x3,y4,'g-', 'LineWidth', 2.00)
grid

axis square
% axis([1 xl 0 8])
%

title('Dependence of standard error on I_1 -- 2.36 bits', 'FontSize', fs1)
%
xlabel('D_1', 'FontSize', fs1)
% ylabel('\nu', 'FontSize', fs1)
ylabel('Standard error, \sigma', 'FontSize', fs1)

xtxt = 2
ytxt = 18
%
text (xtxt, ytxt, ['  NW info = ' num2str(smlde2nw,4), ',  D(1) = ' num2str(dr(1),4),...
',  \sigma = ' num2str(smlde2scnw,4)], 'FontSize', fs2)
text (xtxt, ytxt-5, [' Skew = ' num2str(skewnw,4),...
',  Curt = ' num2str(kurtnw,4)], 'FontSize', fs2)


figure
xl=10000000
%
semilogx (dr,kurt,'g-o', dr,skew,'b-o', dr,sig,'m-o', x3,y2a,'m-',...
x3,y3,'b-', x3,y4,'g-', 'LineWidth', 2.00)

grid

axis square
axis([1 xl -6 8])
%

title('Fig. 9', 'FontSize', fs1)
%title('Dependence of \sigma, \eta, \psi on I_0 -- 2.36 bits', 'FontSize', fs1)
%
xlabel('Derivative |D_1|', 'FontSize', fs1)
% ylabel('\nu', 'FontSize', fs1)
ylabel('Standard error, \sigma; Coefficients \eta, \psi', 'FontSize', fs1)

xtxt = 2
ytxt = 6.5
%
% text (xtxt, ytxt, [' NW info = ' num2str(smlde2nw,4), ',  D_1 = ' num2str(dr(1),4),...
% ',  \sigma = ' num2str(smlde2scnw,4)], 'FontSize', fs2)
% text (xtxt, ytxt-1, [' Skew (\eta_3) = ' num2str(skewnw,4),...
% ',  Kurt (\eta_4) = ' num2str(kurtnw,4)], 'FontSize', fs2)