-
Notifications
You must be signed in to change notification settings - Fork 5
/
allen9.m
155 lines (141 loc) · 4.58 KB
/
allen9.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
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