function allen9(sels) % ALLEN9(sels) % % Magnitude prediction from scale 5 ONLY % as a function of the TriNet magnitude, % the number of stations, and the number of stations versus % TriNet magnitude % % INPUT: % % sels 0 No distance selection % 1 Distance selection to within 100 km % % Last modified by fjsimons-at-alum.mit.edu, 18.07.2006 defval('sels',0) stuff=load('/home/fjsimons/EALARMS/wtALL_diagnostics'); % Collect the data mom=stuff(:,4); dist=stuff(:,5); evnu(:,1)=stuff(:,6); skal(:,1)=stuff(:,7); skal(:,2)=stuff(:,8); skal(:,3)=stuff(:,9); skal(:,4)=stuff(:,10); skal(:,5)=stuff(:,11); % Protect identical moments belonging to individual events % Best regression line is when you DO NOT do this % so use the not-done regression but plot them for yes-done events % We thus define the best regression possible, but show the results % on individual events, not averaging events with the same magnitude mom=mom+evnu/1e10; % Thresholding division thresh=stuff(:,3); difer(diff(thresh)) disp(sprintf('Thresholding divisor %2.2f',thresh(1,1))) if sels==1 % Take out those outside of 100 km selkt=dist<100; mom=mom(selkt); skal=skal(selkt,:); end clf ah=krijetem(subnum(3,1)); slev=20; % Just calculate the magnitude errors at SCALE 5 index=5; % Get the statistics on the logarithmically scaled data [s,g]=row2stats(mom,log10(skal(:,index))); % Use the average regression line on the individual magnitudes instead Gmeanl=[0.9638 -0.4546] % Magnitude break point disp(sprintf('Magnitude breakpoint %4.2f',s(slev))) yfit=polyval(Gmeanl,s(1:slev)); % Figure out the error in the determination % Actual seismic magnitude is in s(1:slev) % Predicted moment is, from the averaged coefficients and the % regression line; predict ALL with the LOW line predmoml=1/Gmeanl(1)*g.mean-Gmeanl(2)/Gmeanl(1); % Make sure this corresponds with what's plotted... disp(sprintf('m_l = %4.2f logC_%i %+4.1f',... 1/Gmeanl(1),index,... -Gmeanl(2)/Gmeanl(1))) % Plot the error versus the TriNet magnitude % Predicting the LOW magnitudes with the LOW regression errll=s(1:slev)-predmoml(1:slev); % Predicting the HIGH magnitudes with the LOW regression errlh=s(slev+1:end)-predmoml(slev+1:end); % Now the fit to the data for the high magnitudes Gmeanh=[0.6851 0.8416] disp(sprintf('m_h= %4.2f logC_%i %+4.1f',... 1/Gmeanh(1),index,... -Gmeanh(2)/Gmeanh(1))) yfit=polyval(Gmeanh,s(slev+1:end)); % Figure out the error in the determination % Predict ALL with the high line predmomh=1/Gmeanh(1)*g.mean-Gmeanh(2)/Gmeanh(1); % Plot the error; % Predicting the LOW magnitudes with the HIGH regression errhl=s(1:slev)-predmomh(1:slev); % Predicting the HIGH magnitudes with the HIGH regression errhh=s(slev+1:end)-predmomh(slev+1:end); % Now the average of both regression lines avhl=mean([predmoml(:) predmomh(:)],2)'; % ERROR VERSUS MAGNITUDE axes(ah(1)) ell=plot(s(1:slev),errll,'s'); hold on [perb,k,l]=errorxy(s,s-avhl,[],abs([errhl errhh]-(s-avhl))); hold on elh=plot(s(slev+1:end),errlh,'s'); ehl=plot(s(1:slev),errhl,'^'); ehh=plot(s(slev+1:end),errhh,'^'); % Labels, etc xl=xlabel('TriNet magnitude'); yl=ylabel('mag prediction error'); grid on pav=plot(s,s-avhl,'o','Markerf','k','MarkerS',3,'MarkerE','k'); % This error IS symmetric, it's about their average hold off % Cosmetics hl=2; set(ehl,'MarkerF',grey,'MarkerS',hl,'MarkerE',grey) set(elh,'MarkerF',grey,'MarkerS',hl,'MarkerE',grey) set(ehh,'MarkerF','k','MarkerS',hl,'MarkerE','k') set(ell,'MarkerF','k','MarkerS',hl,'MarkerE','k') %delete([ehl elh ehh ell]) %delete(pav) delete(perb) set(ah(1),'ylim',[-1.3 1.3],'ytick',[-1:1],'xlim',[2.8 7.4]) % NUMBER OF STATIONS VERSUS TRINET MAGNITUDE axes(ah(2)) nsvm=plot(s,g.nonnans,'d','MarkerS',hl+1,'MarkerE','k','MarkerF','k'); set(ah(2),'ylim',[-5 70],'ytick',... [0:20:70],'xlim',[2.8 7.4]) grid on xl=xlabel('TriNet magnitude'); yl=ylabel('# of reporting stations'); % ERROR VERSUS NUMBER OF REPORTING STATIONS axes(ah(3)) pav=plot(g.nonnans,abs(s-avhl),'o','Markerf','k','MarkerS',3,'MarkerE','k'); % Labels, etc xl=xlabel('number of reporting stations'); yl=ylabel('|mag prediction error|'); grid on % Cosmetics hl=2; set(ah(3),'ylim',[-0.3 1.3],'ytick',[-1:1],'xlim',[0 max(g.nonnans)+1]) fig2print(gcf,'portrait') longticks(ah,2) % Labels [bh1,th1]=label(ah(1),'lr',12); % Move this label about a bit set(bh1,'verti',get(bh1,'verti')+[repmat(0.075,4,1) repmat(0,4,2)]) moveh(th1,0.075) [bh2,th2]=label(ah(2:3),'ur',12,1); set(bh2(1),'verti',get(bh2(1),'verti')+[repmat(0.075,4,1) repmat(0,4,2)]) set(bh2(2),'verti',get(bh2(2),'verti')+... [repmat(1.06,4,1) repmat(0,4,2)]) moveh(th2(1),0.075) moveh(th2(2),1.06) figdisp