% earthquakes_Japan.m, modeled on Bostonquakes.m
% Modeled after LMFig010101poissr.m & earthquakes_Haiti.m
% Data from USGS/NEIC (PDE) 1973 - 2010 12 31 
% http://earthquake.usgs.gov/earthquakes/eqarchives/epic/epic_rect.php
% Search all years
% Japan 9.0 epicenter 38.22 to 142.369 (searched 2.5 n, s e & w
% All quakes greater than 2.0
% Using example 4.3.5 as a model
% Example 4.3.5 Memphis earthquakes and the Poisson distribution
% Pp 305-306 in
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written Fall 2010 by Eugene Gallagher for EEOS601 Revised: 1/9/11
% Eugene.Gallagher@umb.edu 
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% Calculate the exact probability that 2 or more Richter 5 earthquakes will
% strike next year using the Poisson distribution and compare that with the
% normal approximation based on the central limit theorem.
japaneq; % calls the data file.
Richter=4;   % Could be set to any value greater than 3
year=38.2; % Jan 1973 to 3/16/2011
% k=find(DATA>=Richter & DATA<9);
% s=sprintf('Japanese Earthquakes >= Richter Scale %3.1f & R < 9',Richter);
k=find(DATA>=Richter);
s=sprintf('Japanese Earthquakes >= Richter Scale %3.1f',Richter);
edges=Richter:.5:ceil(max(DATA(k)));
N=histc(DATA(k),edges);bar(edges,N/year,'histc');
set(get(gca,'Children'),'FaceColor',[.8 .8 1])
title(s,'FontSize',20);
xlabel('Richter Scale','FontSize',16)
ylabel('Frequency per year','FontSize',16)
ax1=gca;
figure(gcf);pause
i=find(~~N);  % find only the positive values for N.
Y=log(N(i)/year);
X=[ones(length(i),1) edges(i)'];
B=X\Y;
Yest=X*B;
% Now use Matlab's regress to produce the same regression coefficients:
[b,bint,r,rint,stats] = regress(Y,X);
% Note for other problems, more significant figures will be needed for
% reporting the regression coefficients and 95%% CI's
fprintf('\nThe Y intercept = %6.3f with 95%% CI = [%6.3f %6.3f]\n',...
    b(1),bint(1,1),bint(1,2));
fprintf('The slope is %6.3f with 95%% CI = [%6.3f %6.3f]\n',...
    b(2),bint(2,1),bint(2,2));
fprintf(...
'R^2 = %4.1f%%; F = %5.1f with p-value = %7.2g; RMSE = s^2 = %6.3g\n',...
    stats(1)*100,stats(2),stats(3),stats(4))
P=2*(1-tcdf(sqrt(stats(2)),df));
fprintf('P(t >= %4.1f |Ho: slope = 0 & %3.0f df) = %7.2g\n',...
    sqrt(stats(2)),df,P)
semilogy(edges,N/year,'or',edges(i),exp(Yest),'-m');
axis([Richter-0.1    9.1 0.001 109])
xlabel('Richter Scale','FontSize',16)
ylabel('Frequency per year','FontSize',16);
s2=sprintf('N=%8.2f e^{%6.3f R}',exp(B(1)),B(2));
text(6,10,s2,'FontSize',20);
title(s,'FontSize',20)
figure(gcf);pause

% Redo the analysis for earthquakes greater than 4 with log_10 scaling to
% see whether the slope is -1 to match the figure 
% in the text.
Y=log10(N(i)/year);
X=[ones(length(i),1) edges(i)'];
B=X\Y;
Yest=X*B;

% Now use Matlab's regress to produce the same regression coefficients:
[b,bint,r,rint,stats] = regress(Y,X);
% Note for other problems, more significant figures will be needed for
% reporting the regression coefficients and 95%% CI's
fprintf('\nThe Y intercept = %6.3f with 95%% CI = [%6.3f %6.3f]\n',...
    b(1),bint(1,1),bint(1,2));
fprintf('The slope is %6.3f with 95%% CI = [%6.3f %6.3f]\n',...
    b(2),bint(2,1),bint(2,2));
fprintf(...
'R^2 = %4.1f%%; F = %5.1f with p-value = %7.2g; RMSE = s^2 = %6.3g\n',...
    stats(1)*100,stats(2),stats(3),stats(4))
P=2*(1-tcdf(sqrt(stats(2)),df));
fprintf('P(t >= %4.1f |Ho: slope = 0 & %3.0f df) = %7.2g\n',...
    sqrt(stats(2)),df,P)

semilogy(edges,N/year,'or',edges(i),10.^Yest,'-m');
axis([Richter-0.1    9.1 0.009 109])
xlabel('Richter Scale','FontSize',16)
ylabel('Frequency per year','FontSize',16);
s2=sprintf('N=%8.2f 10^{%6.3f R}',10.^B(1),B(2));
text(6,10,s2,'FontSize',20);
title(s,'FontSize',20)
figure(gcf);pause
hold off

semilogy(edges,N/year,'om',edges(i),10.^Yest,'m','LineWidth',3);
set(gca,'ytick',[0.0001 0.001 0.01 0.1 1 10 100])
axis([Richter-0.1    9.1 0.00007 109])
xlabel('Richter Scale','FontSize',16)
ylabel('Frequency per year','FontSize',16);
s2=sprintf('N=%8.2f 10^{%6.3f R}',10.^B(1),B(2));
text(6,10,s2,'FontSize',20);
title(s,'FontSize',20)
hold on
% Continue the line from R 7.5 to  R 9
X=[1 7.5;1 9];
Yest=10.^(X*B);
semilogy(X(:,2),Yest,'m','LineWidth',3);
grid
figure(gcf);pause


Richter=4; 
year=38; % Jan 1973 to 3/16/2011
k=find(DATA>=Richter);
s=sprintf('Japanese Earthquakes >= Richter Scale %3.1f',Richter);
edges=Richter:.5:ceil(max(DATA(k)));
N=histc(DATA(k),edges);
i=find(~~N);  % find only the positive values for N.
Y=log10(N(i)/year);
X=[ones(length(i),1) edges(i)'];
B=X\Y;
Yest2=X*B;
semilogy(edges,N/year,'or',edges(i),10.^Yest2,'r','LineWidth',3);
set(gca,'ytick',[0.0001 0.001 0.01 0.1 1 10 100])
axis([Richter-0.1    9.1 0.00007 109])
figure(gcf)
pause
hold off

