507 lines
16 KiB
Matlab
Executable File
507 lines
16 KiB
Matlab
Executable File
% [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err,BIAS,SD]=
|
|
% Nonpar_tr(t,M,iop,Mmin,Mmax,Nsynth)
|
|
%
|
|
% BASED ON MAGNITUDE SAMPLE DATA M DETERMINES THE ROUND-OFF INTERVAL LENGTH
|
|
% OF THE MAGNITUDE DATA - eps, THE SMOOTHING FACTOR - h, CONSTRUCTS
|
|
% THE BACKGROUND SAMPLE - xx, CALCULATES THE WEIGHTING FACTORS - amb, AND
|
|
% THE END-POINT OF MAGNITUDE DISTRIBUTION Mmax FOR A USE OF THE NONPARAMETRIC
|
|
% ADAPTATIVE KERNEL ESTIMATORS OF MAGNITUDE DISTRIBUTION UNDER THE
|
|
% ASSUMPTION OF THE EXISTENCE OF THE UPPER LIMIT OF MAGNITUDE DISTRIBUTION.
|
|
%
|
|
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UPPER-BOUNDED
|
|
% NON-PARAMETRIC HAZARD ESTIMATION MODE !!
|
|
%
|
|
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
|
|
%
|
|
% DESCRIPTION: The kernel estimator approach is a model-free alternative
|
|
% to estimating the magnitude distribution functions. The smoothing factor
|
|
% h, is estimated using the least-squares cross-validation for the Gaussian
|
|
% kernel function. The final form of the kernel is the adaptive kernel.
|
|
% In order to avoid repetitions, which cannot appear in a sample when the
|
|
% kernel estimators are used, the magnitude sample data are randomized
|
|
% within the magnitude round-off interval. The round-off interval length -
|
|
% eps is the least non-zero difference between sample data or 0.1 is the
|
|
% least difference if greater than 0.1. The randomization is done
|
|
% assuming exponential distribution of m in [m0-eps/2, m0+eps/2], where m0
|
|
% is the sample data point and eps is the length of roud-off inteval. The
|
|
% shape parameter of the exponential distribution is estimated from the whole
|
|
% data sample assuming the exponential distribution. The background sample
|
|
% - xx comprises the randomized values of magnitude doubled symmetrically
|
|
% with respect to the value Mmin-eps/2: length(xx)=2*length(M). Weigthing
|
|
% factors row vector for the adaptive kernel is of the same size as xx.
|
|
% The mean activity rate, lamb, is the number of events >=Mmin into the
|
|
% length of the period in which they occurred.
|
|
% The upper limit of the distribution Mmax is evaluated using
|
|
% the Kijko-Sellevol generic formula. If convergence is not reached the
|
|
% Whitlock @ Robson simplified formula is used:
|
|
% Mmaxest= 2(max obs M) - (second max obs M)).
|
|
%
|
|
% See: the references below for a more comprehensive description.
|
|
%
|
|
% This is a beta version of the program. Further developments are foreseen.
|
|
%
|
|
% REFERENCES:
|
|
%Silverman B.W. (1986) Density Estimation for Statistics and Data Analysis,
|
|
% Chapman and Hall, London
|
|
%Kijko A., Lasocki S., Graham G. (2001) Pure appl. geophys. 158, 1655-1665
|
|
%Lasocki S., Orlecka-Sikora B. (2008) Tectonophysics 456, 28-37
|
|
%Kijko, A., and M.A. Sellevoll (1989) Bull. Seismol. Soc. Am. 79, 3,645-654
|
|
%Lasocki, S., Urban, P. (2011) Acta Geophys 59, 659-673,
|
|
% doi: 10.2478/s11600-010-0049-y
|
|
%
|
|
% INPUT:
|
|
% t - vector of earthquake occurrence times
|
|
% M - vector of earthquake magnitudes (sample data)
|
|
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
|
|
% iop=2 - 'year'
|
|
% Mmin - lower bound of the distribution - catalog completeness level
|
|
% Mmax - upper limit of Magnitude Distribution. Can be set by user, or
|
|
% estimate within the program - it then should be set as Mmax=[].
|
|
%
|
|
% OUTPUT
|
|
% lamb_all - mean activity rate for all events
|
|
% lamb - mean activity rate for events >= Mmin
|
|
% lamb_err - error paramter on the number of events >=Mmin. lamb_err=0
|
|
% for 50 or more events >=Mmin and the parameter estimation is
|
|
% continued, lamb_err=1 otherwise, all output paramters except
|
|
% lamb_all and lamb are set to zero and the function execution is
|
|
% terminated.
|
|
% unit - string with name of time unit used ('year' or 'month' or 'day').
|
|
% eps - length of round-off interval of magnitudes.
|
|
% ierr - h-convergence indicator. ierr=0 if the estimation procedure of
|
|
% the optimal smoothing factor has converged (a zero of the h functional
|
|
% has been found), ierr=1 when multiple zeros of h functional were
|
|
% encountered - the largest h is accepted, ierr = 2 when h functional did
|
|
% not zeroe - the approximate h value is taken.
|
|
% h - kernel smoothing factor.
|
|
% xx - the background sample for the nonparametric estimators of magnitude
|
|
% distribution
|
|
% ambd - the weigthing factors for the adaptive kernel
|
|
% Mmax - upper limit of magnitude distribution
|
|
% err - error parameter on Mmax estimation, err=0 - convergence, err=1 -
|
|
% no converegence of Kijko-Sellevol estimator, Robinson @ Whitlock
|
|
% method used.
|
|
% BIAS - Mmax estimation Bias (Lasocki and Urban, 2011)
|
|
% SD - Mmax standard deviation (Lasocki ands Urban, 2011)
|
|
%
|
|
% 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.
|
|
%
|
|
|
|
function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err,BIAS,SD]=...
|
|
Nonpar_tr_Ob(t,M,iop,Mmin,Mmax,Nsynth)
|
|
if nargin==5;Nsynth=[];end % K08DEC2019
|
|
if isempty(t) || numel(t)<3 isempty(M(M>=Mmin)) %K03OCT
|
|
t=[1 2];M=[1 2]; end %K30SEP
|
|
|
|
lamb_err=0;
|
|
n=length(M);
|
|
t1=t(1);
|
|
%%% %%%%%%%%%%%%%MICHAL
|
|
xx=NaN;
|
|
ambd=NaN;
|
|
%%% %%%%%%%%%%%%%MICHAL
|
|
for i=1:n
|
|
if M(i)>=Mmin; break; end
|
|
t1=t(i+1);
|
|
end
|
|
t2=t(n);
|
|
for i=n:1
|
|
if M(i)>=Mmin; break; end
|
|
t2=t(i-1);
|
|
end
|
|
nn=0;
|
|
for i=1:n
|
|
if M(i)>=Mmin
|
|
nn=nn+1;
|
|
end
|
|
end
|
|
% SL 03MAR2015 ----------------------------------
|
|
[NM,unit]=time_diff(t(1),t(n),iop);
|
|
lamb_all=n/NM;
|
|
[NM,unit]=time_diff(t1,t2,iop);
|
|
lamb=nn/NM;
|
|
% SL 03MAR2015 ----------------------------------
|
|
if nn<50
|
|
eps=0;ierr=0;h=0;Mmax=0;err=0;
|
|
lamb_err=1;
|
|
BIAS=NaN;SD=NaN; %%% K 08NOV2019
|
|
return;
|
|
end
|
|
|
|
eps=magn_accur(M);
|
|
n=0;
|
|
for i=1:length(M)
|
|
if M(i)>=Mmin;
|
|
n=n+1;
|
|
x(n)=M(i);
|
|
end
|
|
end
|
|
x=sort(x)';
|
|
beta=1/(mean(x)-Mmin+eps/2);
|
|
[xx]=korekta(x,Mmin,eps,beta);
|
|
xx=sort(xx);
|
|
clear x;
|
|
xx = podwajanie(xx,Mmin-eps/2);
|
|
[h,ierr]=hopt(xx);
|
|
[ambd]=scaling(xx,h);
|
|
|
|
if isempty(Mmax) %K30AUG2019 - Allow for manually set Mmax
|
|
err_Mmax=1;[Mmax,err]=Mmaxest(xx,h,Mmin-eps/2); % err_Mmax added 04DEC2019
|
|
else
|
|
err_Mmax=0;err=0; %K30AUG2019
|
|
end
|
|
|
|
|
|
|
|
% Estimation of Mmax Bias %%% K 04DEC2019
|
|
% (Lasocki and Urban, 2011, doi:10.2478/s11600-010-0049-y)
|
|
if isempty(Nsynth)==0 && err_Mmax==1 % set number of synthetic datasets, e.g. 10000
|
|
[BIAS,SD]=Mmax_Bias_GR(t,M,Mmin,Mmax,err,h,xx,ambd,Nsynth);
|
|
elseif isempty(Nsynth)==0 && err_Mmax==0;
|
|
warning('Mmax must be empty for BIAS calculation');BIAS=[];SD=[];
|
|
else; BIAS=0;SD=0;
|
|
end
|
|
|
|
end
|
|
|
|
function [NM,unit]=time_diff(t1,t2,iop) % SL 03MAR2015 ----------------------------------
|
|
|
|
% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
|
|
%
|
|
% t1 - start time (in MATLAB numerical format)
|
|
% t2 - end time (in MATLAB numerical format) t2>=t1
|
|
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
|
|
% iop=2 - 'year'
|
|
%
|
|
% NM - number of time units from t1 to t2
|
|
% unit - string with name of time unit used ('year' or 'month' or 'day').
|
|
|
|
if iop==0
|
|
NM=(t2-t1);
|
|
unit='day';
|
|
elseif iop==1
|
|
V1=datevec(t1);
|
|
V2=datevec(t2);
|
|
NM=V2(3)/eomday(V2(1),V2(2))+V2(2)+12-V1(2)-V1(3)/eomday(V1(1),V1(2))...
|
|
+(V2(1)-V1(1)-1)*12;
|
|
unit='month';
|
|
else
|
|
V1=datevec(t1);
|
|
V2=datevec(t2);
|
|
NM2=V2(3);
|
|
if V2(2)>1
|
|
for k=1:V2(2)-1
|
|
NM2=NM2+eomday(V2(1),k);
|
|
end
|
|
end
|
|
day2=365; if eomday(V2(1),2)==29; day2=366; end;
|
|
NM2=NM2/day2;
|
|
NM1=V1(3);
|
|
if V1(2)>1
|
|
for k=1:V1(2)-1
|
|
NM1=NM1+eomday(V1(1),k);
|
|
end
|
|
end
|
|
day1=365; if eomday(V1(1),2)==29; day1=366; end;
|
|
NM1=(day1-NM1)/day1;
|
|
NM=NM2+NM1+V2(1)-V1(1)-1;
|
|
unit='year';
|
|
end
|
|
end
|
|
|
|
function [m_corr]=korekta(m,Mmin,eps,beta)
|
|
|
|
% RANDOMIZATION OF MAGNITUDE WITHIN THE ACCURACY INTERVAL
|
|
%
|
|
% m - input vector of magnitudes
|
|
% Mmin - catalog completeness level
|
|
% eps - accuracy of magnitude
|
|
% beta - the parameter of the unbounded exponential distribution
|
|
%
|
|
% m_corr - vector of randomized magnitudes
|
|
%
|
|
F1=1-exp(-beta*(m-Mmin-0.5*eps));
|
|
F2=1-exp(-beta*(m-Mmin+0.5*eps));
|
|
u=rand(size(m));
|
|
w=u.*(F2-F1)+F1;
|
|
m_corr=Mmin-log(1-w)./beta;
|
|
end
|
|
|
|
function x2 = podwajanie(x,x0)
|
|
|
|
% DOUBLES THE SAMPLE
|
|
|
|
% If the sample x(i) is is truncated from the left hand side and belongs
|
|
% to the interval [x0,inf) or it is truncated from the right hand side and
|
|
% belongs to the interval (-inf,x0]
|
|
% then the doubled sample is [-x(i)+2x0,x(i)]
|
|
% x - is the column data vector
|
|
% x2 - is the column vector of data doubled and sorted in the ascending
|
|
% order
|
|
|
|
x2=[-x+2*x0
|
|
x];
|
|
x2=sort(x2);
|
|
end
|
|
|
|
function [h,ierr]=hopt(x)
|
|
|
|
%Estimation of the optimal smoothing factor by means of the least squares
|
|
%method
|
|
% x - column data vector
|
|
% The result is an optimal smoothing factor
|
|
% ierr=0 - convergence, ierr=1 - multiple h, ierr=2 - approximate h is used
|
|
% The function calls the procedure FZERO for the function 'funct'
|
|
% NEW VERSION 2 - without a square matrix. Also equipped with extra zeros
|
|
% search
|
|
|
|
% MODIFIED JUNE 2014
|
|
|
|
ierr=0;
|
|
n=length(x);
|
|
x=sort(x);
|
|
interval=[0.000001 2*std(x)/n^0.2];
|
|
x1=funct(interval(1),x);
|
|
x2=funct(interval(2),x);
|
|
if x1*x2<0
|
|
fun = @(t) funct(t,x); % for octave
|
|
x0 =interval; % for octave
|
|
[hh(1),fval,exitflag] = fzero(fun,x0); % for octave
|
|
|
|
% Extra zeros search
|
|
jj=1;
|
|
for kk=2:7
|
|
interval(1)=1.1*hh(jj);
|
|
interval(2)=interval(1)+(kk-1)*hh(jj);
|
|
x1=funct(interval(1),x);
|
|
x2=funct(interval(2),x);
|
|
if x1*x2<0
|
|
jj=jj+1;
|
|
fun = @(t) funct(t,x); % for octave
|
|
x0 =interval; % for octave
|
|
[hh(jj),fval,exitflag] = fzero(fun,x0); % for octave
|
|
end
|
|
end
|
|
if jj>1;ierr=1;end
|
|
h=max(hh);
|
|
|
|
if exitflag==1;return;end
|
|
|
|
end
|
|
h=0.891836*(mean(x)-x(1))/(n^0.2);
|
|
ierr=2;
|
|
end
|
|
|
|
function [fct]=funct(t,x)
|
|
p2=1.41421356;
|
|
n=length(x);
|
|
yy=zeros(size(x));
|
|
for i=1:n,
|
|
xij=(x-x(i)).^2/t^2;
|
|
y=exp(-xij/4).*((xij/2-1)/p2)-2*exp(-xij/2).*(xij-1);
|
|
yy(i)=sum(y);
|
|
end;
|
|
fct=sum(yy)-2*n;
|
|
clear xij y yy;
|
|
end
|
|
|
|
|
|
function [ambd]=scaling(x,h)
|
|
|
|
% EVALUATES A VECTOR OF SCALING FACTORS FOR THE NONPARAMETRIC ADAPTATIVE
|
|
% ESTIMATION
|
|
|
|
% x - the n dimensional column vector of data values sorted in the ascending
|
|
% order
|
|
% h - the optimal smoothing factor
|
|
% ambd - the resultant n dimensional row vector of local scaling factors
|
|
|
|
n=length(x);
|
|
c=sqrt(2*pi);
|
|
gau=zeros(1,n);
|
|
for i=1:n,
|
|
gau(i)=sum(exp(-0.5*((x(i)-x)/h).^2))/c/n/h;
|
|
end
|
|
g=exp(mean(log(gau)));
|
|
ambd=sqrt(g./gau);
|
|
end
|
|
|
|
function [eps]=magn_accur(M)
|
|
x=sort(M);
|
|
d=x(2:length(x))-x(1:length(x)-1);
|
|
eps=min(d(d>0));
|
|
if eps>0.1; eps=0.1;end
|
|
end
|
|
|
|
function [Mmax,ierr]=Mmaxest(x,h,Mmin)
|
|
|
|
% ESTIMATION OF UPPER BOUND USING NONPARAMETRIC DISTRIBUTION FUNCTIONS
|
|
% x - row vector of magnitudes (basic sample).
|
|
% h - optimal smoothing factor
|
|
% Mmax - upper bound
|
|
% ierr=0 if basic procedure converges, ierr=1 when Robsen & Whitlock Mmas
|
|
% estimation
|
|
|
|
% Uses function 'dystryb'
|
|
|
|
n=length(x);
|
|
ierr=1;
|
|
x=sort(x);
|
|
Mmax1=x(n);
|
|
for i=1:50,
|
|
d=normcdf((Mmin-x)./h);
|
|
mian=sum(normcdf((Mmax1-x)./h)-d);
|
|
Mmax=x(n)+moja_calka(@dystryb,x(1),Mmax1,0.00001,h,mian,x,d);
|
|
if abs(Mmax-Mmax1)<0.01
|
|
ierr=0;break;
|
|
end
|
|
Mmax1=Mmax;
|
|
end
|
|
if (ierr==1 || Mmax>9)
|
|
Mmax=2*x(n)-x(n-1);
|
|
ierr=1;
|
|
end
|
|
end
|
|
|
|
function [y]=dystryb(z,h,mian,x,d)
|
|
n=length(x);
|
|
m=length(z);
|
|
for i=1:m,
|
|
t=(z(i)-x)./h;
|
|
t=normcdf(t);
|
|
yy=sum(t-d);
|
|
y(i)=(yy/mian)^n;
|
|
end
|
|
end
|
|
|
|
function [calka,ier]=moja_calka(funfc,a,b,eps,varargin)
|
|
|
|
% Integration by means of 16th poit Gauss method. Adopted from CERNLIBRARY
|
|
|
|
% funfc - string with the name of function to be integrated
|
|
% a,b - integration limits
|
|
% eps - accurracy
|
|
% varargin - other parameters of function to be integrated
|
|
% calka - integral
|
|
% ier=0 - convergence, ier=1 - no conbergence
|
|
|
|
persistent W X CONST
|
|
W=[0.101228536290376 0.222381034453374 0.313706645877887 ...
|
|
0.362683783378362 0.027152459411754 0.062253523938648 ...
|
|
0.095158511682493 0.124628971255534 0.149595988816577 ...
|
|
0.169156519395003 0.182603415044924 0.189450610455069];
|
|
X=[0.960289856497536 0.796666477413627 0.525532409916329 ...
|
|
0.183434642495650 0.989400934991650 0.944575023073233 ...
|
|
0.865631202387832 0.755404408355003 0.617876244402644 ...
|
|
0.458016777657227 0.281603550779259 0.095012509837637];
|
|
CONST=1E-12;
|
|
delta=CONST*abs(a-b);
|
|
calka=0.;
|
|
aa=a;
|
|
y=b-aa;
|
|
ier=0;
|
|
while abs(y)>delta,
|
|
bb=aa+y;
|
|
c1=0.5*(aa+bb);
|
|
c2=c1-aa;
|
|
s8=0.;
|
|
s16=0.;
|
|
for i=1:4,
|
|
u=X(i)*c2;
|
|
s8=s8+W(i)*(feval(funfc,c1+u,varargin{:})+feval(funfc,c1-u,varargin{:}));
|
|
end
|
|
for i=5:12,
|
|
u=X(i)*c2;
|
|
s16=s16+W(i)*(feval(funfc,c1+u,varargin{:})+feval(funfc,c1-u,varargin{:}));
|
|
end
|
|
s8=s8*c2;
|
|
s16=s16*c2;
|
|
if abs(s16-s8)>eps*(1+abs(s16))
|
|
y=0.5*y;
|
|
calka=0.;
|
|
ier=1;
|
|
else
|
|
calka=calka+s16;
|
|
aa=bb;
|
|
y=b-aa;
|
|
ier=0;
|
|
end
|
|
end
|
|
end
|
|
|
|
% --------------------- Mmax BIAS estimation routine ---------------------- K 08NOV2019
|
|
|
|
function [BIAS,SD]=Mmax_Bias_GR(t,m,Mc,Mmax1,err,h,xx,ambd,synth)
|
|
|
|
s1=sort(unique(m));s2=s1(2:length(s1))-s1(1:length(s1)-1);EPS=min(s2);
|
|
|
|
if err~=0
|
|
warning('process did not converge!!')
|
|
end
|
|
|
|
MAXm=max(m);N=numel(m(m>=Mc));DeltaM=MAXm-Mc; %beta=b*log(10);
|
|
[mag,PDF_NPT,CDF_NPT]=dist_NPT(Mc-EPS,Mmax1+EPS,0.01,Mc,EPS,h,xx,ambd,Mmax1);
|
|
|
|
for j=1:synth %set number of synthetic datasets, default is 10000
|
|
% % CDF:
|
|
|
|
% j
|
|
|
|
% linear interpolation to assign magnitude values from a uniform distribution sample
|
|
iM=rand(1,N);M1=interp1q(CDF_NPT,mag,iM');
|
|
br(j)=1/(log(10)*(mean(M1)-min(M1)+EPS/2));DM=range(M1);
|
|
Mmax=max(M1);
|
|
|
|
% Iteration Process to estimate b and Mmax
|
|
b1=10;best=[1.0 10.0];i=1;
|
|
while min(abs(diff(best)))>0.00001
|
|
w=exp(b1*(Mmax-Mc));E1=expint(N/(w-1));E2=expint(N*w/(w-1));
|
|
%E=expint(w);
|
|
Mme=Mmax+(E1-E2)/(b1*exp(-N/(w-1)))+(Mc)*exp(-N); %Mme=round(Mme/EPS)*EPS;
|
|
|
|
if isnan(Mme)
|
|
KM=sort(unique(M1),'descend');
|
|
Mme=2*KM(1)-KM(2);
|
|
end
|
|
|
|
fun=@(bb) 1/bb+(Mme-Mc)/(1-exp(bb*(Mme-Mc)))-mean(M1)+Mc-EPS/2; %consider th5 last 0.05 term
|
|
b1=fzero(fun,1);best(i)=b1;i=i+1;
|
|
|
|
if i==50
|
|
warning('process did not converge!!');break
|
|
end
|
|
|
|
end
|
|
|
|
be(j)=b1/log(10);
|
|
Me(j)=Mme;dm(j)=DM;Mm(j)=Mmax;
|
|
|
|
end
|
|
|
|
BIAS=mean(MAXm-Me);
|
|
SD=std(MAXm-Me);
|
|
|
|
%b-mean(be) %check b-value difference
|
|
%histogram(be)
|
|
|
|
% MAXm: maximum magnitude in the real catalog
|
|
% Mmax: maximum magnitudes observed in the synthetic catalogs (rounded)
|
|
% Me: maximum magnitude estimates for the synthetic catalogs
|
|
% Mmax1: maximum magnitude estimated by GRT
|
|
end
|
|
|
|
|
|
|
|
|