Added MagnitudeConversion scripts
This commit is contained in:
parent
e866b3d977
commit
944922ba75
7
src/MagnitudeConversion/LICENCE.txt
Normal file
7
src/MagnitudeConversion/LICENCE.txt
Normal file
@ -0,0 +1,7 @@
|
||||
This application is a part of the IS-EPOS e-PLATFORM.
|
||||
|
||||
This is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.
|
||||
186
src/MagnitudeConversion/LinReg.m
Normal file
186
src/MagnitudeConversion/LinReg.m
Normal file
@ -0,0 +1,186 @@
|
||||
% function [msyn,a,b,sa,sb,cc,rms,N]=LinReg(m1,m2,opt,n)
|
||||
% -----------------------------------------------------------------------------------------------
|
||||
% The program calculates the linear regression parameters (a and b)
|
||||
% and their uncertainties by applying different regression techniques.
|
||||
% -----
|
||||
% The program calls the following functions, appended in the scipt:
|
||||
% function [a,b,sa,sb,rms,N]=GOR(m1,m2,n) - General Orthogonal Regression
|
||||
% function [a b sa sb rms]=OLS(m1,m2) - Ordinary Least Squares
|
||||
% function [a b sa sb rms]=ILS(m1,m2) - Inverted Least Squares
|
||||
% -----------------------------------------------------------------------------------------------
|
||||
% COMPLIED BY:
|
||||
% Leptokaropoulos Konstantinos, June 2014, (kleptoka@igf.edu.pl) within IS-EPOS Project,
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
% DESCRIPTION
|
||||
% The program is used to establish relationships between different magnitude scales.
|
||||
% These relationships can be then used for converting diverse magnitude scales into a
|
||||
% common one, such that a homogeneous magnitude is applied for further analyses.
|
||||
% The suggested method is GOR, because the ordinary least-squares assumes that
|
||||
% there are no uncertainties in the values of the independent variable, a condition that
|
||||
% is never truth when dealing with earthquake magnitudes. This may introduce
|
||||
% systematic errors in magnitude conversion, apparent catalog incompleteness, and
|
||||
% significant bias in the estimates of the b-value. According to GOR technique, the
|
||||
% projection of the independent variable is done along a weighted orthogonal distance
|
||||
% from the linear fitting curve. Although the error varriance ratio of the two variables
|
||||
% is introduced in calculations, Castellaro et al. (2006) showed that even if the applied
|
||||
% values differ from the real ones, GOR method still performs better than the OLS.
|
||||
% Despite the GOR superiority, OLS and ILS are also provided as supplementary
|
||||
% options in the present program, for comparison.
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
% INPUT PARAMETERS:
|
||||
% m1 - variable 1 (e.g. magnitude 1) [VECTOR]
|
||||
% m2 - variable 2 (e.g. magnitude 2) [VECTOR]
|
||||
% opt - regression method selection (1,2 or 3): [SCALAR] (Default=1)
|
||||
% 1 - General Orthogonal Regression - GOR function
|
||||
% 2 - Ordinary Least Squares - OLS function
|
||||
% 3 - Inverted Least Squares - ILS function
|
||||
% n - error variance ratio (var(m2)/var(m1)) [SCALAR] (FOR GOR only, Default=1)
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
% OUTPUT:
|
||||
% msyn - Converted Magnitude ([VECTOR]
|
||||
% a - intercept [SCALAR]
|
||||
% b - slope [SCALAR]
|
||||
% sa - 95% confidence interval for intercept [SCALAR]
|
||||
% sb - 95% confidence interval for slope [SCALAR]
|
||||
% cc - Pearson correlation coefficient between vectors m1 and m2 [SCALAR-optional]
|
||||
% rms - root mean square error [SCALAR-optional]
|
||||
% N - number of data pairs m1,m2 [SCALAR-optional]
|
||||
% ---------------------------------------------------------------------------------------------
|
||||
% REFERENCES
|
||||
% **** (for General Orthogonal Regression) ****
|
||||
% Fuller,W. A. (1987). Measurement Error Models,Wiley, New York, 440 pp.
|
||||
% Castellaro, S., F. Mulargia, and Y. Y. Kagan (2006). Regression problems for magnitudes,
|
||||
% Geophys. J. Int. 165, 913–930.
|
||||
% Lolli, B., and P. Gasperini (2012). A comparison among general orthogonal regression methods
|
||||
% applied to earthquake magnitude conversions, Geophys. J. Int. 190, 1135–1151.
|
||||
% Wason, H. R., R. Das, and M. L. Sharma (2012).Magnitude conversion problem using general
|
||||
% orthogonal regression, Geophys. J. Int. 190, 1091–1096.
|
||||
% ----------------------------------------------------------------------------------------------
|
||||
% LICENSE
|
||||
% This file is a part of the IS-EPOS e-PLATFORM.
|
||||
%
|
||||
% This is free software: you can redistribute it and/or modify it under
|
||||
% the terms of the GNU General Public License as published by the Free
|
||||
% Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% This program is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
|
||||
function [msyn,a,b,sa,sb,cc,rms,N]=LinReg(m1,m2,opt,n)
|
||||
|
||||
%%%%% Remove NaN values - K03JUN2015 %%%%%
|
||||
mm1=m1(isnan(m1)==0 & isnan(m2)==0);
|
||||
mm2=m2(isnan(m1)==0 & isnan(m2)==0);
|
||||
m1=mm1;m2=mm2;
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
|
||||
if nargin==2, opt=1; n=1;
|
||||
elseif nargin==3 , n=1;end
|
||||
|
||||
% CC
|
||||
cc=corrcoef(m1,m2);
|
||||
cc=cc(2);
|
||||
N=length(m1);
|
||||
|
||||
if opt==1
|
||||
[a,b,sa,sb,rms]=GOR(m1,m2,n);
|
||||
elseif opt==2
|
||||
[a,b,sa,sb,rms]=OLS(m1,m2);
|
||||
elseif opt==3
|
||||
[a,b,sa,sb,rms]=ILS(m1,m2);
|
||||
end
|
||||
msyn=m1*b+a;
|
||||
end
|
||||
|
||||
% ********************************************************
|
||||
% ------------------------- GOR --------------------------------
|
||||
function [a,b,sa,sb,rms,N]=GOR(m1,m2,n)
|
||||
|
||||
if nargin==2 || isempty(n),n=1;end
|
||||
|
||||
%GOR parameters
|
||||
f=1; %alternative f=0
|
||||
cova=cov(m1,m2,f);
|
||||
mx=var(m1,f);
|
||||
my=var(m2,f);
|
||||
mxy=cova(1,2);
|
||||
|
||||
b=(my-n*mx+sqrt((my-n*mx)^2+4*n*mxy^2))/(2*mxy);
|
||||
a=mean(m2)-b*mean(m1);
|
||||
|
||||
% BASIC STATISTICS
|
||||
msyn=b*m1+a;
|
||||
rms=sqrt(sum((m2-msyn).^2)/(length(m2)));
|
||||
N=length(m1);
|
||||
|
||||
% Confidence Intervals of a and b-values
|
||||
d=(my-n*mx);
|
||||
|
||||
sx=(sqrt(d^2+4*n*mxy^2)-d)/(2*n);
|
||||
sy=(my+n*mx-sqrt(d^2+4*n*mxy^2))/(2*n);
|
||||
|
||||
s=((N-1)*(n+b^2)*sy)/(N-2);
|
||||
|
||||
Vb=(mx*s-(b*sy)^2)/((N-1)*sx^2);
|
||||
Va=s/N+Vb*mean(m1)^2;
|
||||
|
||||
sb=sqrt(Vb)*tinv(0.975,N-2);
|
||||
sa=sqrt(Va)*tinv(0.975,N-2);
|
||||
end
|
||||
|
||||
% ----------------------------- OLS ----------------------------
|
||||
function [a b sa sb rms]=OLS(m1,m2)
|
||||
f=1;N=length(m1);
|
||||
cova=cov(m1,m2,f);
|
||||
mx=var(m1,f);
|
||||
my=var(m2,f);
|
||||
mxy=cova(1,2);
|
||||
|
||||
p=polyfit(m1,m2,1);
|
||||
|
||||
vb=my/sum((m1-mean(m1)).^2);
|
||||
va=my*sum(m1.^2)/(N*sum((m1-mean(m1)).^2));
|
||||
sb=tinv(0.975,N-2)*sqrt(vb);
|
||||
sa=tinv(0.975,N-2)*sqrt(va);
|
||||
b=p(1);a=p(2);
|
||||
|
||||
msyn=b*m1+a;
|
||||
rms=sqrt(sum((m2-msyn).^2)/(length(m2)));
|
||||
end
|
||||
|
||||
% ----------------------------- ILS ----------------------------
|
||||
function [a b sa sb rms]=ILS(m1,m2)
|
||||
f=1;N=length(m1);
|
||||
cova=cov(m1,m2,f);
|
||||
mx=var(m2,f);
|
||||
my=var(m1,f);
|
||||
mxy=cova(1,2);
|
||||
|
||||
p=polyfit(m2,m1,1);
|
||||
q=[1/p(1) -p(2)/p(1)];
|
||||
|
||||
|
||||
vb=my/sum((m2-mean(m2)).^2);
|
||||
va=my*sum(m2.^2)/(N*sum((m2-mean(m2)).^2));
|
||||
sb=tinv(0.975,N-2)*sqrt(vb);
|
||||
sa=tinv(0.975,N-2)*sqrt(va);
|
||||
|
||||
sq1=abs(sb*q(1)/p(1));
|
||||
sq2=abs(sb*q(2)/p(2));
|
||||
b=q(1);a=q(2);sb=sq1;sa=sq2;
|
||||
|
||||
msyn=b*m1+a;
|
||||
rms=sqrt(sum((m2-msyn).^2)/(length(m2)));
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
186
src/MagnitudeConversion/LinReg_O.m
Normal file
186
src/MagnitudeConversion/LinReg_O.m
Normal file
@ -0,0 +1,186 @@
|
||||
% function [msyn,a,b,sa,sb,cc,rms,N]=LinReg_O(m1,m2,opt,n) - Octave Compatible Version
|
||||
% -----------------------------------------------------------------------------------------------
|
||||
% The program calculates the linear regression parameters (a and b)
|
||||
% and their uncertainties by applying different regression techniques.
|
||||
% -----
|
||||
% The program calls the following functions, appended in the scipt:
|
||||
% function [a,b,sa,sb,rms,N]=GOR(m1,m2,n) - General Orthogonal Regression
|
||||
% function [a b sa sb rms]=OLS(m1,m2) - Ordinary Least Squares
|
||||
% function [a b sa sb rms]=ILS(m1,m2) - Inverted Least Squares
|
||||
% -----------------------------------------------------------------------------------------------
|
||||
% COMPLIED BY:
|
||||
% Leptokaropoulos Konstantinos, June 2014, (kleptoka@igf.edu.pl) within IS-EPOS Project,
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
% DESCRIPTION
|
||||
% The program is used to establish relationships between different magnitude scales.
|
||||
% These relationships can be then used for converting diverse magnitude scales into a
|
||||
% common one, such that a homogeneous magnitude is applied for further analyses.
|
||||
% The suggested method is GOR, because the ordinary least-squares assumes that
|
||||
% there are no uncertainties in the values of the independent variable, a condition that
|
||||
% is never truth when dealing with earthquake magnitudes. This may introduce
|
||||
% systematic errors in magnitude conversion, apparent catalog incompleteness, and
|
||||
% significant bias in the estimates of the b-value. According to GOR technique, the
|
||||
% projection of the independent variable is done along a weighted orthogonal distance
|
||||
% from the linear fitting curve. Although the error varriance ratio of the two variables
|
||||
% is introduced in calculations, Castellaro et al. (2006) showed that even if the applied
|
||||
% values differ from the real ones, GOR method still performs better than the OLS.
|
||||
% Despite the GOR superiority, OLS and ILS are also provided as supplementary
|
||||
% options in the present program, for comparison.
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
% INPUT PARAMETERS:
|
||||
% m1 - variable 1 (e.g. magnitude 1) [VECTOR]
|
||||
% m2 - variable 2 (e.g. magnitude 2) [VECTOR]
|
||||
% opt - regression method selection (1,2 or 3): [SCALAR] (Default=1)
|
||||
% 1 - General Orthogonal Regression - GOR function
|
||||
% 2 - Ordinary Least Squares - OLS function
|
||||
% 3 - Inverted Least Squares - ILS function
|
||||
% n - error variance ratio (var(m2)/var(m1)) [SCALAR] (FOR GOR only, Default=1)
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
% OUTPUT:
|
||||
% msyn - Converted Magnitude ([VECTOR]
|
||||
% a - intercept [SCALAR]
|
||||
% b - slope [SCALAR]
|
||||
% sa - 95% confidence interval for intercept [SCALAR]
|
||||
% sb - 95% confidence interval for slope [SCALAR]
|
||||
% cc - Pearson correlation coefficient between vectors m1 and m2 [SCALAR-optional]
|
||||
% rms - root mean square error [SCALAR-optional]
|
||||
% N - number of data pairs m1,m2 [SCALAR-optional]
|
||||
% ---------------------------------------------------------------------------------------------
|
||||
% REFERENCES
|
||||
% **** (for General Orthogonal Regression) ****
|
||||
% Fuller,W. A. (1987). Measurement Error Models,Wiley, New York, 440 pp.
|
||||
% Castellaro, S., F. Mulargia, and Y. Y. Kagan (2006). Regression problems for magnitudes,
|
||||
% Geophys. J. Int. 165, 913–930.
|
||||
% Lolli, B., and P. Gasperini (2012). A comparison among general orthogonal regression methods
|
||||
% applied to earthquake magnitude conversions, Geophys. J. Int. 190, 1135–1151.
|
||||
% Wason, H. R., R. Das, and M. L. Sharma (2012).Magnitude conversion problem using general
|
||||
% orthogonal regression, Geophys. J. Int. 190, 1091–1096.
|
||||
% ----------------------------------------------------------------------------------------------
|
||||
% LICENSE
|
||||
% This file is a part of the IS-EPOS e-PLATFORM.
|
||||
%
|
||||
% This is free software: you can redistribute it and/or modify it under
|
||||
% the terms of the GNU General Public License as published by the Free
|
||||
% Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% This program is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
% ------------------------------------------------------------------------------------------------
|
||||
|
||||
function [msyn,a,b,sa,sb,cc,rms,N]=LinReg_O(m1,m2,opt,n)
|
||||
|
||||
%%%%% Remove NaN values - K03JUN2015 %%%%%
|
||||
mm1=m1(isnan(m1)==0 & isnan(m2)==0);
|
||||
mm2=m2(isnan(m1)==0 & isnan(m2)==0);
|
||||
m1=mm1;m2=mm2;
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
if nargin==2, opt=1; n=1;
|
||||
elseif nargin==3 , n=1;end
|
||||
|
||||
% CC
|
||||
cc=corr(m1,m2);
|
||||
N=length(m1);
|
||||
|
||||
if opt==1
|
||||
[a,b,sa,sb,rms]=GOR(m1,m2,n);
|
||||
elseif opt==2
|
||||
[a,b,sa,sb,rms]=OLS(m1,m2);
|
||||
elseif opt==3
|
||||
[a,b,sa,sb,rms]=ILS(m1,m2);
|
||||
end
|
||||
msyn=m1*b+a;
|
||||
end
|
||||
|
||||
% ********************************************************
|
||||
% ------------------------- GOR --------------------------------
|
||||
function [a,b,sa,sb,rms,N]=GOR(m1,m2,n)
|
||||
|
||||
if nargin==2 || isempty(n),n=1;end
|
||||
|
||||
%GOR parameters
|
||||
f=1; %alternative f=0
|
||||
mx=var(m1,f);
|
||||
my=var(m2,f);
|
||||
mxy=cov(m1,m2,f);
|
||||
%mxy=mxy(2); Enable for Matlab version
|
||||
|
||||
b=(my-n*mx+sqrt((my-n*mx)^2+4*n*mxy^2))/(2*mxy);
|
||||
a=mean(m2)-b*mean(m1);
|
||||
|
||||
% BASIC STATISTICS
|
||||
msyn=b*m1+a;
|
||||
rms=sqrt(sum((m2-msyn).^2)/(length(m2)));
|
||||
N=length(m1);
|
||||
|
||||
% Confidence Intervals of a and b-values
|
||||
d=(my-n*mx);
|
||||
|
||||
sx=(sqrt(d^2+4*n*mxy^2)-d)/(2*n);
|
||||
sy=(my+n*mx-sqrt(d^2+4*n*mxy^2))/(2*n);
|
||||
|
||||
s=((N-1)*(n+b^2)*sy)/(N-2);
|
||||
|
||||
Vb=(mx*s-(b*sy)^2)/((N-1)*sx^2);
|
||||
Va=s/N+Vb*mean(m1)^2;
|
||||
|
||||
sb=sqrt(Vb)*tinv(0.975,N-2);
|
||||
sa=sqrt(Va)*tinv(0.975,N-2);
|
||||
end
|
||||
|
||||
% ----------------------------- OLS ----------------------------
|
||||
function [a b sa sb rms]=OLS(m1,m2)
|
||||
f=1;N=length(m1);
|
||||
mxy=cov(m1,m2,f);
|
||||
mx=var(m1,f);
|
||||
my=var(m2,f);
|
||||
|
||||
p=polyfit(m1,m2,1);
|
||||
|
||||
vb=my/sum((m1-mean(m1)).^2);
|
||||
va=my*sum(m1.^2)/(N*sum((m1-mean(m1)).^2));
|
||||
sb=tinv(0.975,N-2)*sqrt(vb);
|
||||
sa=tinv(0.975,N-2)*sqrt(va);
|
||||
b=p(1);a=p(2);
|
||||
|
||||
msyn=b*m1+a;
|
||||
rms=sqrt(sum((m2-msyn).^2)/(length(m2)));
|
||||
end
|
||||
|
||||
% ----------------------------- ILS ----------------------------
|
||||
function [a b sa sb rms]=ILS(m1,m2)
|
||||
f=1;N=length(m1);
|
||||
mxy=cov(m1,m2,f);
|
||||
mx=var(m1,f);
|
||||
my=var(m2,f);
|
||||
|
||||
p=polyfit(m2,m1,1);
|
||||
q=[1/p(1) -p(2)/p(1)];
|
||||
|
||||
|
||||
vb=my/sum((m2-mean(m2)).^2);
|
||||
va=my*sum(m2.^2)/(N*sum((m2-mean(m2)).^2));
|
||||
sb=tinv(0.975,N-2)*sqrt(vb);
|
||||
sa=tinv(0.975,N-2)*sqrt(va);
|
||||
|
||||
sq1=abs(sb*q(1)/p(1));
|
||||
sq2=abs(sb*q(2)/p(2));
|
||||
b=q(1);a=q(2);sb=sq1;sa=sq2;
|
||||
|
||||
msyn=b*m1+a;
|
||||
rms=sqrt(sum((m2-msyn).^2)/(length(m2)));
|
||||
end
|
||||
|
||||
% Testing plotting
|
||||
%plot(c.Ml,c.Mw,'ko')
|
||||
%hold on
|
||||
%plot([min(c.Ml) max(c.Ml)],[min(c.Ml)*b+a max(c.Ml)*b+a],'r--')
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user