tcs-ah/sera-applications
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update distribution functions for Mmax bias

Browse Source
This commit is contained in:
Konstantinos Leptokaropoulos
2020-06-03 08:36:41 +02:00
parent e7412096b1
commit b4e572d7f2
50 changed files with 7680 additions and 7366 deletions
Download Patch File Download Diff File Expand all files Collapse all files

99
SHAPE_Package/SHAPE_ver1.0/SSH/Nonpar_tr_O.m
View File

@@ -1,5 +1,5 @@
% [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err]=
% Nonpar_tr(t,M,iop,Mmin)
% [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
@@ -55,6 +55,8 @@
% 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
@@ -79,6 +81,8 @@
% 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.
@@ -94,9 +98,9 @@
% GNU General Public License for more details.
%
function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err]=...
Nonpar_tr_O(t,M,iop,Mmin,Mmax)
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
@@ -131,6 +135,7 @@ lamb=nn/NM;
if nn<50
eps=0;ierr=0;h=0;Mmax=0;err=0;
lamb_err=1;
BIAS=NaN;SD=NaN; %%% K 08NOV2019
return;
end
@@ -152,15 +157,21 @@ xx = podwajanie(xx,Mmin-eps/2);
[ambd]=scaling(xx,h);
if isempty(Mmax) %K30AUG2019 - Allow for manually set Mmax
[Mmax,err]=Mmaxest(xx,h,Mmin-eps/2);
err_Mmax=1;[Mmax,err]=Mmaxest(xx,h,Mmin-eps/2); % err_Mmax added 04DEC2019
else
err=0; %K30AUG2019
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
% enai=dlmread('paraT.txt'); %for fixed xx,ambd to test in different platforms
% [ambd]=enai(:,1);
% xx=enai(:,2)';
% [h,ierr]=hopt(xx);
% [Mmax,err]=Mmaxest(xx,h,Mmin-eps/2);
end
@@ -429,3 +440,67 @@ while abs(y)>delta,
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

92
SHAPE_Package/SHAPE_ver1.0/SSH/TruncGR_O.m
View File

@@ -1,5 +1,5 @@
%
% [lamb_all,lamb,lmab_err,unit,eps,b,Mmax,err]=TruncGR(t,M,iop,Mmin)
%[lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_Ob(t,M,iop,Mmin,Mmax,Nsynth)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE,
@@ -38,6 +38,8 @@
% iop=2 - 'year'
% Mmin - catalog completeness level. Must be determined externally.
% Can take any value from [min(M), max(M)].
% Mmax - upper limit of Magnitude Distribution. Can be set by user, or
% estimate within the program - it then should be set as Mmax=[].
%
% OUTPUT:
%
@@ -55,6 +57,8 @@
% 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.
@@ -70,7 +74,8 @@
% GNU General Public License for more details.
%
function [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err]=TruncGR_O(t,M,iop,Mmin,Mmax)
function [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_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
@@ -103,6 +108,7 @@ lamb=nn/NM;
if nn<15
eps=0;b=0;Mmax=0;err=0;
lamb_err=1;
BIAS=NaN;SD=NaN; %%% K 08NOV2019
return;
end
@@ -121,14 +127,16 @@ Max_obs=max(xx);
beta0=0;
Mmax1=Max_obs;
if isempty(Mmax)==0 %%% K 28JUL2015
err_Mmax=0; %%% K 04DEC2019
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax); %%% K 28JUL2015
x0 = 1; %[0.05,4.0]; %%% K 28JUL2015 - See exception line 153
x0 = 1; %[0.05,4.0]; %%% K 28JUL2015 - See exception line 155
beta = fzero(fun,x0); %%% K 28JUL2015
err=0; %%% K 28JUL2015
else %%% K 28JUL2015 - line 148
else %%% K 28JUL2015 - line 150
err_Mmax=1; %%% K 04DEC2019
for i=1:50,
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax1);
x0 =1; %[0.05,4.0]; %%% K29JUL2015 - See exception line 153
x0 =1; %[0.05,4.0]; %%% K29JUL2015 - See exception line 155
beta = fzero(fun,x0);
Mmax=Max_obs+moja_calka('f_podc',Mmin,Max_obs,1e-5,nn,beta,Mmin-eps/2,Mmax1);
if ((abs(Mmax-Mmax1)<0.01)&&(abs(beta-beta0)<0.0001))
@@ -152,8 +160,22 @@ clear xx
% Exception for v-value
if b<0.05 || b>6.0; error('Unacceptable b-value, abort and select different dataset');end
beta;
% 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,b,err,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
@@ -303,3 +325,63 @@ d=x(2:length(x))-x(1:length(x)-1);
eps=min(d(d>0));
if eps>0.1; eps=0.1;end
end
% --------------------- Mmax BIAS estimation routine ---------------------- K 08NOV2019
function [BIAS,SD]=Mmax_Bias_GR(t,m,Mc,Mmax1,b,err,synth)
if err~=0
warning('process did not converge!!'),pause
end
MAXm=max(m);beta=b*log(10);N=numel(m(m>=Mc));DeltaM=MAXm-Mc;
for j=1:synth %set number of synthetic datasets, default is 10000
% % CDF:
M=Mc:0.0001:MAXm;upt=1-exp(-beta*(M-Mc));
dwt=1-exp(-beta*(MAXm-Mc));F=upt./dwt; % j
% linear interpolation to assign magnitude values from a uniform distribution sample
iM=rand(1,N);M1=interp1q(F',M',iM');
br(j)=1/(log(10)*(mean(M1)-min(M1)));DM=range(M1);
Mmax=max(M1);
% Iteration Process to estimate b and Mmax
b1=1;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; %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

99
SHAPE_Package/SHAPE_ver2.0/SSH/Nonpar_tr_O.m
View File

@@ -1,5 +1,5 @@
% [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err]=
% Nonpar_tr(t,M,iop,Mmin)
% [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
@@ -55,6 +55,8 @@
% 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
@@ -79,6 +81,8 @@
% 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.
@@ -94,9 +98,9 @@
% GNU General Public License for more details.
%
function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err]=...
Nonpar_tr_O(t,M,iop,Mmin,Mmax)
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
@@ -131,6 +135,7 @@ lamb=nn/NM;
if nn<50
eps=0;ierr=0;h=0;Mmax=0;err=0;
lamb_err=1;
BIAS=NaN;SD=NaN; %%% K 08NOV2019
return;
end
@@ -152,15 +157,21 @@ xx = podwajanie(xx,Mmin-eps/2);
[ambd]=scaling(xx,h);
if isempty(Mmax) %K30AUG2019 - Allow for manually set Mmax
[Mmax,err]=Mmaxest(xx,h,Mmin-eps/2);
err_Mmax=1;[Mmax,err]=Mmaxest(xx,h,Mmin-eps/2); % err_Mmax added 04DEC2019
else
err=0; %K30AUG2019
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
% enai=dlmread('paraT.txt'); %for fixed xx,ambd to test in different platforms
% [ambd]=enai(:,1);
% xx=enai(:,2)';
% [h,ierr]=hopt(xx);
% [Mmax,err]=Mmaxest(xx,h,Mmin-eps/2);
end
@@ -429,3 +440,67 @@ while abs(y)>delta,
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

92
SHAPE_Package/SHAPE_ver2.0/SSH/TruncGR_O.m
View File

@@ -1,5 +1,5 @@
%
% [lamb_all,lamb,lmab_err,unit,eps,b,Mmax,err]=TruncGR(t,M,iop,Mmin)
%[lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_Ob(t,M,iop,Mmin,Mmax,Nsynth)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE,
@@ -38,6 +38,8 @@
% iop=2 - 'year'
% Mmin - catalog completeness level. Must be determined externally.
% Can take any value from [min(M), max(M)].
% Mmax - upper limit of Magnitude Distribution. Can be set by user, or
% estimate within the program - it then should be set as Mmax=[].
%
% OUTPUT:
%
@@ -55,6 +57,8 @@
% 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.
@@ -70,7 +74,8 @@
% GNU General Public License for more details.
%
function [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err]=TruncGR_O(t,M,iop,Mmin,Mmax)
function [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_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
@@ -103,6 +108,7 @@ lamb=nn/NM;
if nn<15
eps=0;b=0;Mmax=0;err=0;
lamb_err=1;
BIAS=NaN;SD=NaN; %%% K 08NOV2019
return;
end
@@ -121,14 +127,16 @@ Max_obs=max(xx);
beta0=0;
Mmax1=Max_obs;
if isempty(Mmax)==0 %%% K 28JUL2015
err_Mmax=0; %%% K 04DEC2019
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax); %%% K 28JUL2015
x0 = 1; %[0.05,4.0]; %%% K 28JUL2015 - See exception line 153
x0 = 1; %[0.05,4.0]; %%% K 28JUL2015 - See exception line 155
beta = fzero(fun,x0); %%% K 28JUL2015
err=0; %%% K 28JUL2015
else %%% K 28JUL2015 - line 148
else %%% K 28JUL2015 - line 150
err_Mmax=1; %%% K 04DEC2019
for i=1:50,
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax1);
x0 =1; %[0.05,4.0]; %%% K29JUL2015 - See exception line 153
x0 =1; %[0.05,4.0]; %%% K29JUL2015 - See exception line 155
beta = fzero(fun,x0);
Mmax=Max_obs+moja_calka('f_podc',Mmin,Max_obs,1e-5,nn,beta,Mmin-eps/2,Mmax1);
if ((abs(Mmax-Mmax1)<0.01)&&(abs(beta-beta0)<0.0001))
@@ -152,8 +160,22 @@ clear xx
% Exception for v-value
if b<0.05 || b>6.0; error('Unacceptable b-value, abort and select different dataset');end
beta;
% 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,b,err,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
@@ -303,3 +325,63 @@ d=x(2:length(x))-x(1:length(x)-1);
eps=min(d(d>0));
if eps>0.1; eps=0.1;end
end
% --------------------- Mmax BIAS estimation routine ---------------------- K 08NOV2019
function [BIAS,SD]=Mmax_Bias_GR(t,m,Mc,Mmax1,b,err,synth)
if err~=0
warning('process did not converge!!'),pause
end
MAXm=max(m);beta=b*log(10);N=numel(m(m>=Mc));DeltaM=MAXm-Mc;
for j=1:synth %set number of synthetic datasets, default is 10000
% % CDF:
M=Mc:0.0001:MAXm;upt=1-exp(-beta*(M-Mc));
dwt=1-exp(-beta*(MAXm-Mc));F=upt./dwt; % j
% linear interpolation to assign magnitude values from a uniform distribution sample
iM=rand(1,N);M1=interp1q(F',M',iM');
br(j)=1/(log(10)*(mean(M1)-min(M1)));DM=range(M1);
Mmax=max(M1);
% Iteration Process to estimate b and Mmax
b1=1;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; %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