% earthquakes_LA.m, modeled on Bostonquakes.m
% See Also earthquakes_Japan.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 1/1/73 to 12/31/10 38 years
% LA is -118.25, 34.05
% LA 33.05 to 35.05 lat Long -119.25 to -117.25
% All quakes greater than 2.0
% a) copy data from screen to Matlab window, b) save as comma delimited
% text file, c) open file in spreadsheet (excel or quatro pro) d)
% delete unnecessary columns, especially last column of characters that
% fouls up a Matlab load. e) paste Richter magnitude data into Matlab 
% editor and save.
% 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: 3/18/11
% Eugene.Gallagher@umb.edu 
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
earthquakes_LA_38y; % calls the data file.
Richter=4;   % Could be set to any value greater than 2 for LA quakes
year=38;
k=find(DATA(:,7)>=4);
edges=Richter:.5:ceil(max(DATA(k,7)));
N=histc(DATA(k,7),edges);bar(edges,N/year,'histc');
s=sprintf('LA Earthquakes >= Richter Scale %3.1f',Richter);
title(s,'FontSize',20);
set(get(gca,'Children'),'FaceColor',[.8 .8 1]);
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    6.7 0.009 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(4.2296,10,s2,'FontSize',20);
title(s,'FontSize',20)
figure(gcf);pause

% Redo the analysis with for earthquakes greater than 4 with log_10 scaling to
% see whether the slope is -1to 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    6.7 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(4.2296,10,s2,'FontSize',20);
title(s,'FontSize',20)
figure(gcf);pause