tcs-ah/is-epos-scripts
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

magnitude completeness script added

Browse Source
This commit is contained in:
Joanna Kocot
2016-02-11 11:00:32 +01:00
parent 54ddf1b8e5
commit fbba8ae4fd
1 changed files with 231 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

231
src/CompletenessMagnitudeEstimation/MCcalc.m Normal file
View File

@@ -0,0 +1,231 @@
% function [OUTstruct,maxc,gft10,gft5,mgft]=MCcalc(m,t,s,iop)
% ----------------------------------------------------------------------------------------------
% Program that calculates the completess magnitude of a given
% catalog by application of different methods, ASSUMING THAT
% FREQUENCY-MANGITUDE DISTRIBUTION IS EXPONENTIAL
% The program calls the following functions which are appended in the scipt:
% function aa=FM_Param(m,t,iop) - GR-law parameters as a function of mmin
% function [b]=bmle(x,y,ml,dm) - Maximum likelihood b-value estimation
% function Rs=rescalc(No,m,a,b) - Residuals calculation
% ----------------------------------------------------------------------------------------------
% DESCRIPTION:
% The completeness magnitude is estimated by the application of Goodness<EFBFBD>of<EFBFBD>Fit Test (GFT),
% proposed by Wiemer and Wyss (2000). The procedure they followed is that a power law,
% as a function of minimum magnitude, Mi, is fitted for events with M?Mi, by application
% of maximum likelihood estimation (Aki, 1965). The synthetic data, i.e. the distribution of
% magnitudes which represent a perfect fit to the power law, is constructed in this way. Then
% the normalized absolute difference, R, between the cumulative number of observed events
% (No) and the simulated ones (Ns) in each magnitude bin is computed. A model is found at
% an R-value at which a predefined percentage (90% and 95%) of the observed data is modeled
% by a straight line, which means that 90% and 95% of the observed data can be simulated by
% the specific power law. In addition to this approach, synthetic catalogs are constructed following
% the defined G-R, law, and then the number of events for each magnitude bin of the synthetic
% catalogs is compared with the corresponding observed one. In this way a discrete minimum
% residual point is calculated, and considered as the completeness magnitude for a given dataset.
% The Maximum curvature method is also applied for comparison, which corresponds to the
% magnitude bin with the highest frequency of events in the non-cumulative F-M distribution.
% NOTE: Recall that the assumption of self similarity of the magnitude distribution may be invalid
% in some cases. In such cases, alternative (e.g. network based) methods should be applied.
% THE USER SHOULD COMPREHEND THE PROGRAM RESULTS TOGETHER WITH THEIR
% VISUALIZATION IN ORDER TO MAKE HIS/HERS FINAL DECISION.
% ----------------------------------------------------------------------------------------------
% COMPLIED BY:
% Leptokaropoulos Konstantinos, last updated December 2015, (kleptoka@igf.edu.pl)
% within IS-EPOS Project,
% ----------------------------------------------------------------------------------------------
% INPUT:
% m: magnitude of the events in the catalog
% t: occurrence time of the events
% s: number of synthetic catalogs following the same distribution with
% the observed data (for modified goodness of fit test)-default is 500
% iop: unit time selection (default is 'day'):
% 0: day
% 1: month
% 2: year
%--------------------------------------------------------------------------------------------
% OUTPUT
% aa: output (not extracted) is an 9-column array with seismicity parameters as
% a function of minimum magnitude. This array is also saved as
% 'parameters.txt'. Its columns correspond to:
% aa(1): magnitude bin (the step is fixed to 0.1)
% aa(2): number of events in the respective magnitude bin
% aa(3): cumulative number of events in the respective magnitude bin
% aa(4): b-value estimated by Maximum Likelihood (Aki, 1965)
% function bmle, appended
% aa(5): b-value uncertainty (Aki, 1965)
% aa(6): a-value calculated from GR law logN=a-bM, given N,M from catalog
% and b from MLE (aa(4)).
% aa(7): at-value, normalized per unit time (at=a-log10(DT))
% aa(8): Residual between observed data and power law,
% as a function of minimum magnitude
% aa(9): Average residual between synthetic catalogs (s) and power law,
% as a function of minimum magnitude
%
% OUTstruct: Structure containing all the aforementioned parameters in matlab format
% maxc: Mc estimation by Maximum Curvature Method(Wiemer & Wyss, 2000)
% gft10: Mc estimation by Goodness of fit test at 90% confidence bounds (Wiemer & Wyss, 2000)
% gft5: Mc estimation by Goodness of fit test at 95% confidence bounds (Wiemer & Wyss, 2000)
% mgft: Mc estimation by Modified goodness of fit test (Leptokaropoulos et al., 2013)
%----------------------------------------------------------------------------------------------
% REFERENCES
% - Aki, K. (1965), Maximum likelihood estimate of b in the formula log N=a-bM and its
% confidence limits, Bull. Earthq. Res. Inst., Univ. Tokyo, 43, 237-239
% - Wiemer, S. and M. Wyss (2000), Minimum magnitude of completeness in earthquake catalogs:
% Examples from Alaska, the Western United States, and Japan, Bull. Seismol. Soc. Am., 90, 859-869.
% - Leptokaropoulos et al., (2013), A homogeneous earthquake catalogue compilation for western
% Turkey and magnitude of completeness determination, Bull. Seismol. Soc. Am., 103, 5, 2739-2751.
%------------------------------------------------------------------------------------------------
function [OUTstruct,maxc,gft10,gft5,mgft]=MCcalc(m,t,s,iop)
if nargin==2;s=500;iop=0;end
if nargin==3;iop=0;end
% -------------- VALIDATION RULES -------------
if s<10 || s>10000;error('Number of synthetic catalogs should be between 10 and 10000');end
%----------------------------------------------------------
%call FM_Param function
aa=FM_Param(m,t,iop);
no=aa(:,2);
% Nobs-Npl residuals
No=aa(:,3);m=aa(:,1);b=aa(:,4);a=aa(:,6);
n=length(No);ml=min(m);mu=max(m);
% call rescalc function
Rs=rescalc(No,m,a,b);
% Nsynth-Npl residuals
for i=1:s
ra=rand(No(1),1);
nN=No/No(1);
te=histc(ra,flipud(nN));
te(1)=numel(find(ra<min(nN)));
te=flipud(cumsum(te));
Sy=te(2:n);Sy(n)=numel(find(ra<min(nN)));
Rsyn=rescalc(Sy,m,a,b);
R(i,:)=Rsyn;
Rf=mean(R(:,1:n))';
end
Rf(Rf==inf)=100; Rf(Rf>100)=100; % K 17APR2014
maxc=max(m(no==max(no)));
gft10=min(m(Rs<10)); gft5=min(m(Rs<5));
mgft=m(Rf==min(Rf)); rmgft=Rf(Rf==min(Rf));
if gft10==max(m);gft10=NaN;end %K4NOE2014 %K 12NOV 2015- replace [] with NaNs
if gft5==max(m);gft5=NaN;end %K4NOE2014 %K 12NOV 2015- replace [] with NaNs
aa(:,8)=Rs;aa(:,9)=Rf;
[OUTstruct M n N b berr a at Rs Rf]=ar2vec(aa);
end
% function aa=FM_Param(m,t,iop)
% Function that calculates and saves seismicity parameters as a
% function of minimum magnitude for given catalog
% -------------------------------------------------------------------------------------------
% -------------------------------------------------------------------------------------------
function aa=FM_Param(m,t,iop)
if nargin==2;iop=0;end
dt=range(t);
if iop==0 % unit time
dt=dt;un='day';
elseif iop==1;
dt=dt/30.4;un='month';
elseif iop==2
dt=dt/365;un='year';
end
dm=0.1;
m=round(m*10)*dm; %rounds magnitides to 1 decimal digit
mu=max(m);ml=min(m); % maximum catalog magnitude
mb=(ml:dm:mu)'; % magnitude bins
n1=hist(m,min(m):0.1:max(m))'; % K-15MAY2015 number of events belonging to each magnitude bin
nx=length(n1);
cn1=flipud(cumsum(n1(nx:-1:1))); % cumulative number of events in each magnitude bin
for i=1:nx % calculation of b-value as a function of magnitude
b1=bmle(mb(i:nx),n1(i:nx),mb(i),dm);
sb(i)=b1/sqrt(cn1(i)); % b-value error
a(i)=log10(cn1(i))+b1*mb(i); % a-value
b(i)=b1;
end
b=b';sb=sb';a=a';
a1=a-log10(dt); % equivalent a-value for unit-time period
aa=[mb n1 cn1 b sb a a1]; % output parameters
end
%------------------------------------------------------------
%------------------------------------------------------------
function [b]=bmle(x,y,ml,dm)
% Function calcutates the b-value of frequency-magnitude
% distribution of earthquakes following the Maximum
% Likelihood Estimate (Aki, 1965).
% INPUT:
% x: magnitude
% y: number of events with magnitude x
% ml: minimum magnitude (cuttoff)
% dm: the binning width of the catalog (default dm=0.1)
me=sum(x.*y)/sum(y); % arithmetic mean of sample
b=1/((me-ml+dm/2)*log(10));
end
%------------------------------------------------------------
%------------------------------------------------------------
function Rs=rescalc(No,m,a,b)
% Residuals Calculation
% Input:
% No: number of observed data per magnitude bin
% m: magnitude vector
% a: GR law a-value
% b: GR law b-value
% OUTPUT
% Rs: Fitting residuals as a function of minimum magnitude
n=length(No);
for i=1:n
Ns=10.^(a(i)-m*b(i));
Rs(i)=sum(abs(Ns(i:n)-No(i:n)))/sum(No(i:n));
end
Rs=Rs'*100;
end
%------------------------------------------------------------
%------------------------------------------------------------ Added 17DEC2015
function [OUTstruct M n N b berr a at Rs Rf]=ar2vec(aa)
% Separates the output result 'aa' of the 'MCcalc'
% function to 9 vectors and 1 structure, corresponding to:
% M: Magnitude bins, from the minimum to the maximum
% magnitude of the data set.
% n: Number of events in the respective magnitude bin
% N: Cumulative number of events in the respective magnitude bin
% b: b-value of the GR law
% berr: b-value uncertainty
% a: a-value of the GR law
% at: at-value (a-value normalized per unit time)
% Rs: Residual between observed data and power law
% Rf: Average residual between synthetic catalogs (s) and power law
%
%OUTstruct: is a structure containing the aforementioned fields
% -----------------------------------------------------------------
% Input is the 'aa' array as it is derived from 'MCcalc' program
%------------------------------------------------------------------
M=aa(:,1);
n= aa(:,2);
N=aa(:,3);
b= aa(:,4);
berr= aa(:,5);
a =aa(:,6);
at=aa(:,7);
Rs= aa(:,8);
Rf=aa(:,9);
OUTstruct.M=M;OUTstruct.n=n;OUTstruct.N=N;
OUTstruct.b=b;OUTstruct.berr=berr;OUTstruct.a=a;
OUTstruct.at=at;OUTstruct.Rs=Rs;OUTstruct.Rf=Rf;
end