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

stationary seismic hazard analysis scripts added

Browse Source
This commit is contained in:
Joanna Kocot
2016-01-12 10:26:08 +01:00
commit 54ddf1b8e5
21 changed files with 2352 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
.gitignore vendored Normal file
View File

@@ -0,0 +1,4 @@
/*
/*/
!.gitignore
!/src/

94
src/StationarySeismicHazardAnalysis/ExcProbGRT.m Normal file
View File

@@ -0,0 +1,94 @@
% [x,z]=ExcProbGRT(opt,xd,xu,dx,y,Mmin,lamb,eps,b,Mmax)
%
%EVALUATES THE EXCEEDANCE PROBABILITY VALUES USING THE UPPER-BOUNDED G-R
% LED MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the upper-bounded Gutenberg-Richter
% relation leads to the upper truncated exponential distribution to model
% magnitude distribution from and above the catalog completness level
% Mmin. The shape parameter of this distribution, consequently the G-R
% b-value and the end-point of the distriobution Mmax as well as the
% activity rate of M>=Mmin events are calculated at start-up of the
% stationary hazard assessment services in the upper-bounded
% Gutenberg-Richter estimation mode.
%
% The exceedance probability of magnitude M' in the time period of
% length T' is the probability of an earthquake of magnitude M' or greater
% to occur in T'. Depending on the value of the parameter opt the
% exceedance probability values are calculated for a fixed time period T'
% and different magnitude values or for a fixed magnitude M' and different
% time period length values. In either case the independent variable vector
% starts from xd, up to xu with step dx. In either case the result is
% returned in the vector z.
%
%INPUT:
% opt - determines the mode of calculations. opt=0 - fixed time period
% length (y), different magnitude values (x), opt=1 - fixed magnitude
% (y), different time period lengths (x)
% xd - starting value of the changeable independent variable
% xu - ending value of the changeable independent variable
% dx - step change of the changeable independent variable
% y - fixed independent variable value: time period length T' if opt=0,
% magnitude M' if opt=1
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
% Mmax - upper limit of magnitude distribution
%OUTPUT:
% x - vector of changeable independent variable: magnitudes if opt=0,
% time period lengths if opt=1,
% x=(xd:dx:xu)
% z - vector of exceedance probability values of the same length as x
%
% 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, <http://www.gnu.org/licenses/>.
%
function [x,z]=ExcProbGRT(opt,xd,xu,dx,y,Mmin,lamb,eps,b,Mmax)
beta=b*log(10);
if opt==0
if xd<Mmin; xd=Mmin;end
if xu>Mmax; xu=Mmax;end
end
x=(xd:dx:xu)';
if opt==0
z=1-exp(-lamb*y.*(1-Cdfgr(x,beta,Mmin-eps/2,Mmax)));
else
z=1-exp(-lamb*(1-Cdfgr(y,beta,Mmin-eps/2,Mmax)).*x);
end
end
function [y]=Cdfgr(t,beta,Mmin,Mmax)
%CDF of the truncated upper-bounded exponential distribution (truncated G-R
% model
% Mmin - catalog completeness level
% Mmax - upper limit of the distribution
% beta - the distribution parameter
% t - vector of magnitudes (independent variable)
% y - CDF vector
mian=(1-exp(-beta*(Mmax-Mmin)));
y=(1-exp(-beta*(t-Mmin)))/mian;
idx=find(y>1);
y(idx)=ones(size(idx));
end

74
src/StationarySeismicHazardAnalysis/ExcProbGRU.m Normal file
View File

@@ -0,0 +1,74 @@
% [x,z]=ExcProbGRU(opt,xd,xu,dx,y,Mmin,lamb,eps,b)
%
%EVALUATES THE EXCEEDANCE PROBABILITY VALUES USING THE UNLIMITED G-R
% LED MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the unlimited Gutenberg-Richter relation
% leads to the exponential distribution model of magnitude distribution
% from and above the catalog completness level Mmin. The shape parameter of
% this distribution and consequently the G-R b-value are calculated at
% start-up of the stationary hazard assessment services in the
% unlimited Gutenberg-Richter estimation mode.
%
% The exceedance probability of magnitude M' in the time period of
% length T' is the probability of an earthquake of magnitude M' or greater
% to occur in T'. Depending on the value of the parameter opt the
% exceedance probability values are calculated for a fixed time period T'
% and different magnitude values or for a fixed magnitude M' and different
% time period length values. In either case the independent variable vector
% starts from xd, up to xu with step dx. In either case the result is
% returned in the vector z.
%
%INPUT:
% opt - determines the mode of calculations. opt=0 - fixed time period
% length (y), different magnitude values (x), opt=1 - fixed magnitude
% (y), different time period lengths (x)
% xd - starting value of the changeable independent variable
% xu - ending value of the changeable independent variable
% dx - step change of the changeable independent variable
% y - fixed independent variable value: time period length T' if opt=0,
% magnitude M' if opt=1
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
%OUTPUT
% x - vector of changeable independent variable: magnitudes if opt=0,
% time period lengths if opt=1,
% x=(xd:dx:xu)
% z - vector of exceedance probability values of the same length as x
%
% 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, <http://www.gnu.org/licenses/>.
%
function [x,z]=ExcProbGRU(opt,xd,xu,dx,y,Mmin,lamb,eps,b)
beta=b*log(10);
if opt==0
if xd<Mmin; xd=Mmin;end
end
x=(xd:dx:xu)';
if opt==0
z=1-exp(-lamb*y.*exp(-beta*(x-Mmin+eps/2)));
else
z=1-exp(-lamb*exp(-beta*(y-Mmin+eps/2)).*x);
end
end

112
src/StationarySeismicHazardAnalysis/ExcProbNPT.m Normal file
View File

@@ -0,0 +1,112 @@
% [x,z]=ExcProbNPT(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd,Mmax)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES THE
% EXCEEDANCE PROBABILITY VALUES FOR THE UPPER-BOUNDED NONPARAMETRIC
% DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution has a hard end point Mmax from the right hand
% side.The estimation makes use of the previously estimated parameters
% namely the mean activity rate lamb, the length of magnitude round-off
% interval, eps, the smoothing factor, h, the background sample, xx, the
% scaling factors for the background sample, ambd, and the end-point of
% magnitude distribution Mmax. The background sample,xx, comprises the
% randomized values of observed magnitude doubled symmetrically with
% respect to the value Mmin-eps/2.
%
% The exceedance probability of magnitude M' in the time
% period of length T' is the probability of an earthquake of magnitude M'
% or greater to occur in T'.
%
% Depending on the value of the parameter opt the exceedance probability
% values are calculated for a fixed time period T' and different magnitude
% values or for a fixed magnitude M' and different time period length
% values. In either case the independent variable vector starts from
% xd, up to xu with step dx. In either case the result is returned in the
% vector z.
%
% 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
%
% INPUT:
% opt - determines the mode of calculations. opt=0 - fixed time period
% length (y), different magnitude values (x), opt=1 - fixed magnitude
% (y), different time period lengths (x)
% xd - starting value of the changeable independent variable
% xu - ending value of the changeable independent variable
% dx - step change of the changeable independent variable
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
% Mmax - upper limit of magnitude distribution
%
% OUTPUT:
% x - vector of changeable independent variable x=(xd:dx:xu)
% z - vector of exceedance probability values
%
% 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, <http://www.gnu.org/licenses/>.
%
function [x,z]=...
ExcProbNPT(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd,Mmax)
if opt==0
if xd<Mmin; xd=Mmin;end
if xu>Mmax; xu=Mmax;end
end
x=(xd:dx:xu)';
n=length(x);
mian=2*(Dystr_npr(Mmax,xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
if opt==0
for i=1:n
CDF_NPT=2*(Dystr_npr(x(i),xx,ambd,h)...
-Dystr_npr(Mmin-eps/2,xx,ambd,h))./mian;
z(i)=1-exp(-lamb*y.*(1-CDF_NPT));
end
else
CDF_NPT=2*(Dystr_npr(y,xx,ambd,h)...
-Dystr_npr(Mmin-eps/2,xx,ambd,h))./mian;
z=1-exp(-lamb*(1-CDF_NPT).*x);
if y>Mmax;z=zeros(size(x));end
end
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

101
src/StationarySeismicHazardAnalysis/ExcProbNPU.m Normal file
View File

@@ -0,0 +1,101 @@
% [x,z]=ExcProbNPU(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES THE
% EXCEEDANCE PROBABILITY VALUES FOR THE UNBOUNDED NONPARAMETRIC
% DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution is unlimited from the right hand side.
% The estimation makes use of the previously estimated parameters of kernel
% estimation, namely the smoothing factor, the background sample and the
% scaling factors for the background sample. The background sample
% - xx comprises the randomized values of observed magnitude doubled
% symmetrically with respect to the value Mmin-eps/2.
% The exceedance probability of magnitude M' in the time period of length
% T' is the probability of an earthquake of magnitude M' or greater to
% occur in T'.
% Depending on the value of the parameter opt the exceedance probability
% values are calculated for a fixed time period T' and different magnitude
% values or for a fixed magnitude M' and different time period length
% values. In either case the independent variable vector starts from
% xd, up to xu with step dx. In either case the result is returned in the
% vector z.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation fro 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
%
% INPUT:
% opt - determines the mode of calculations. opt=0 - fixed time period
% length (y), different magnitude values (x), opt=1 - fixed magnitude
% (y), different time period lengths (x)
% xd - starting value of the changeable independent variable
% xu - ending value of the changeable independent variable
% dx - step change of the changeable independent variable
% y - fixed independent variable value: time period length T' if opt=0,
% magnitude M' if opt=1
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
%
% OUTPUT:
% x - vector of changeable independent variable: magnitudes if opt=0,
% time period lengths if opt=1,
% x=(xd:dx:xu)
% z - vector of exceedance probability values of the same length as x
%
% 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, <http://www.gnu.org/licenses/>.
%
function [x,z]=ExcProbNPU(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd)
x=(xd:dx:xu)';
n=length(x);
if opt==0
for i=1:n
CDF_NPU=2*(Dystr_npr(x(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
z(i)=1-exp(-lamb*y.*(1-CDF_NPU));
end
else
CDF_NPU=2*(Dystr_npr(y,xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
z=1-exp(-lamb*(1-CDF_NPU).*x);
end
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

55
src/StationarySeismicHazardAnalysis/Max_credM_GRT.m Normal file
View File

@@ -0,0 +1,55 @@
% [T,m]=Max_credM_GRT(Td,Tu,dT,Mmin,lamb,eps,b,Mmax)
%EVALUATES THE MAXIMUM CREDIBLE MAGNITUDE VALUES USING THE UPPER-BOUNDED
% G-R LED MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the upper-bounded Gutenberg-Richter
% relation leads to the upper truncated exponential distribution to model
% magnitude distribution from and above the catalog completness level
% Mmin. The shape parameter of this distribution, consequently the G-R
% b-value and the end-point of the distriobution Mmax as well as the
% activity rate of M>=Mmin events are calculated at start-up of the
% stationary hazard assessment services in the upper-bounded
% Gutenberg-Richter estimation mode.
%
% The maximum credible magnitude values are calculated for periods of
% length starting from Td up to Tu with step dT.
%
% INPUT:
% Td - starting period length for maximum credible magnitude calculations
% Tu - ending period length for maximum credible magnitude calculations
% dT - period length step for maximum credible magnitude calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
% Mmax - upper limit of magnitude distribution
%
% OUTPUT:
% T - vector of independent variable (period lengths) T=(Td:dT:Tu)
% m - vector of maximum credible magnitudes of the same length as T
%
% 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, <http://www.gnu.org/licenses/>.
%
function [T,m]=Max_credM_GRT(Td,Tu,dT,Mmin,lamb,eps,b,Mmax)
T=(Td:dT:Tu)';
beta=b*log(10);
mian=(1-exp(-beta*(Mmax-Mmin+eps/2)));
m=Mmin-eps/2-1/beta*log((1-(1-1./(lamb*T))*mian));
end

59
src/StationarySeismicHazardAnalysis/Max_credM_GRU.m Normal file
View File

@@ -0,0 +1,59 @@
% [T,m]=Max_credM_GRU(Td,Tu,dT,Mmin,lamb,eps,b)
%
%EVALUATES THE MAXIMUM CREDIBLE MAGNITUDE VALUES USING THE UNLIMITED
% G-R LED MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the unlimited Gutenberg-Richter relation
% leads to the exponential distribution model of magnitude distribution
% from and above the catalog completness level Mmin. The shape parameter of
% this distribution and consequently the G-R b-value are calculated at
% start-up of the stationary hazard assessment services in the
% unlimited Gutenberg-Richter estimation mode.
%
% The maximum credible magnitude for the period of length T
% is the magnitude value whose mean return period is T.
%
% The maximum credible magnitude values are calculated for periods of
% length starting from Td up to Tu with step dT.
%
%INPUT:
% Td - starting period length for maximum credible magnitude calculations
% Tu - ending period length for maximum credible magnitude calculations
% dT - period length step for maximum credible magnitude calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
%
%OUTPUT:
% T - vector of independent variable (period lengths) T=(Td:dT:Tu)
% m - vector of maximum credible magnitudes of the same length as T
%
%
% 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, <http://www.gnu.org/licenses/>.
%
function [T,m]=Max_credM_GRU(Td,Tu,dT,Mmin,lamb,eps,b)
T=(Td:dT:Tu)';
beta=b*log(10);
m=Mmin-eps/2+1/beta.*log(lamb*T);
end

94
src/StationarySeismicHazardAnalysis/Max_credM_NPT.m Normal file
View File

@@ -0,0 +1,94 @@
% [T,m]=Max_credM_NPT(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd,Mmax)
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES THE MAXIMUM
% CREDIBLE MAGNITUDE VALUES FOR THE UPPER-BOUNDED NONPARAMETRIC
% DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution has a hard end point Mmax from the right hand
% side.The estimation makes use of the previously estimated parameters
% namely the mean activity rate lamb, the length of magnitude round-off
% interval, eps, the smoothing factor, h, the background sample, xx, the
% scaling factors for the background sample, ambd, and the end-point of
% magnitude distribution Mmax. The background sample,xx, comprises the
% randomized values of observed magnitude doubled symmetrically with
% respect to the value Mmin-eps/2.
%
% The maximum credible magnitude for the period of length T
% is the magnitude value whose mean return period is T.
% The maximum credible magnitude values are calculated for periods of
% length starting from Td up to Tu with step dT.
%
% 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
%
% INPUT:
% Td - starting period length for maximum credible magnitude calculations
% Tu - ending period length for maximum credible magnitude calculations
% dT - period length step for maximum credible magnitude calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
% Mmax - upper limit of magnitude distribution
%
% OUTPUT:
% T - vector of independent variable (period lengths) T=(Td:dT:Tu)
% m - vector of maximum credible magnitudes of the same length as T
%
% 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, <http://www.gnu.org/licenses/>.
%
function [T,m]=Max_credM_NPT(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd,Mmax)
T=(Td:dT:Tu)';
n=length(T);
interval=[Mmin-eps/2 Mmax-0.001];
for i=1:n
m(i)=fzero(@F_maxmagn,interval,[],xx,h,ambd,Mmin-eps/2,Mmax,lamb,T(i));
end
m=m';
end
function [y]=F_maxmagn(t,xx,h,ambd,xmin,Mmax,lamb,D)
mian=2*(Dystr_npr(Mmax,xx,ambd,h)-Dystr_npr(xmin,xx,ambd,h));
CDF_NPT=2*(Dystr_npr(t,xx,ambd,h)-Dystr_npr(xmin,xx,ambd,h))/mian;
y=CDF_NPT-1+1/(lamb*D);
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

94
src/StationarySeismicHazardAnalysis/Max_credM_NPU.m Normal file
View File

@@ -0,0 +1,94 @@
% [T,m]=Max_credM_NPU(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES
% THE MAXIMUM CREDIBLE MAGNITUDE VALUES FOR THE UNBOUNDED
% NONPARAMETRIC DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution is unlimited from the right hand side.
% The estimation makes use of the previously estimated parameters of kernel
% estimation, namely the smoothing factor, the background sample and the
% scaling factors for the background sample. The background sample
% - xx comprises the randomized values of observed magnitude doubled
% symmetrically with respect to the value Mmin-eps/2.
%
% The maximum credible magnitude for the period of length T
% is the magnitude value whose mean return period is T.
% The maximum credible magnitude values are calculated for periods of
% length starting from Td up to Tu with step dT.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation fro 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
%
%INPUT:
% opt - determines the mode of calculations. opt=0 - fixed time period
% length (y), different magnitude values (x), opt=1 - fixed magnitude
% (y), different time period lengths (x)
% xd - starting value of the changeable independent variable
% xu - ending value of the changeable independent variable
% dx - step change of the changeable independent variable
% y - fixed independent variable value: time period length T' if opt=0,
% magnitude M' if opt=1
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
%
%OUTPUT:
% T - vector of independent variable (period lengths) T=(Td:dT:Tu)
% m - vector of maximum credible magnitudes of the same length as T
%
% 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, <http://www.gnu.org/licenses/>.
%
function [T,m]=Max_credM_NPU(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd)
T=(Td:dT:Tu)';
n=length(T);
interval=[Mmin-eps/2 10.0];
for i=1:n
m(i)=fzero(@F_maxmagn_NPU,interval,[],xx,h,ambd,Mmin-eps/2,lamb,T(i));
end
m=m';
end
function [y]=F_maxmagn_NPU(t,xx,h,ambd,xmin,lamb,D)
CDF_NPU=2*(Dystr_npr(t,xx,ambd,h)-Dystr_npr(xmin,xx,ambd,h));
y=CDF_NPU-1+1/(lamb*D);
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

257
src/StationarySeismicHazardAnalysis/Nonpar.m Normal file
View File

@@ -0,0 +1,257 @@
% [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd]=Nonpar(t,M,iop,Mmin)
%
% 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 AND CALCULATES THE WEIGHTING FACTORS - ambd
% FOR A USE OF THE NONPARAMETRIC ADAPTATIVE KERNEL ESTIMATORS OF MAGNITUDE
% DISTRIBUTION.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% 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.
% 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
%
% 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
%
% 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 (the 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
%
% 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, <http://www.gnu.org/licenses/>.
%
function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd]=...
Nonpar(t,M,iop,Mmin)
lamb_err=0;
n=length(M);
t1=t(1);
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
if iop==0
lamb_all=n/round(t(n)-t(1));
lamb=nn/round(t2-t1);
unit='day';
elseif iop==1
lamb_all=30*n/(t(n)-t(1));
lamb=30*nn/(t2-t1);
unit='month';
else
lamb_all=365*n/(t(n)-t(1));
lamb=365*nn/(t2-t1);
unit='year';
end
if nn<50
eps=0;ierr=0;h=0;
lamb_err=1;
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 = doubling(xx,Mmin-eps/2);
[h,ierr]=hopt(xx);
[ambd]=scaling(xx,h);
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 = doubling(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
[hh(1),fval,exitflag]=fzero(@funct,interval,[],x);
% 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;
[hh(jj),fval,exitflag]=fzero(@funct,interval,[],x);
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

371
src/StationarySeismicHazardAnalysis/Nonpar_tr.m Normal file
View File

@@ -0,0 +1,371 @@
% [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err]=
% Nonpar(t,M,iop,Mmin)
%
% 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.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% 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
%
% 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.
%
% 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, <http://www.gnu.org/licenses/>.
%
function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd,Mmax,err]=...
Nonpar_tr(t,M,iop,Mmin)
lamb_err=0;
n=length(M);
t1=t(1);
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
if iop==0
lamb_all=n/round(t(n)-t(1));
lamb=nn/round(t2-t1);
unit='day';
elseif iop==1
lamb_all=30*n/(t(n)-t(1));
lamb=30*nn/(t2-t1);
unit='month';
else
lamb_all=365*n/(t(n)-t(1));
lamb=365*nn/(t2-t1);
unit='year';
end
if nn<50
eps=0;ierr=0;h=0;Mmax=0;err=0;
lamb_err=1;
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 = doubling(xx,Mmin-eps/2);
[h,ierr]=hopt(xx);
[ambd]=scaling(xx,h);
[Mmax,err]=Mmaxest(xx,h,Mmin-eps/2);
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 = doubling(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
[hh(1),fval,exitflag]=fzero(@funct,interval,[],x);
% 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;
[hh(jj),fval,exitflag]=fzero(@funct,interval,[],x);
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

79
src/StationarySeismicHazardAnalysis/Ret_periodGRT.m Normal file
View File

@@ -0,0 +1,79 @@
% [m,T]=Ret_periodGRT(Md,Mu,dM,Mmin,lamb,eps,b,Mmax)
%
% EVALUATES THE MEAN RETURN PERIOD VALUES USING THE UPPER-BOUNDED G-R LED
% MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the upper-bounded Gutenberg-Richter
% relation leads to the upper truncated exponential distribution to model
% magnitude distribution from and above the catalog completness level
% Mmin. The shape parameter of this distribution, consequently the G-R
% b-value and the end-point of the distriobution Mmax as well as the
% activity rate of M>=Mmin events are calculated at start-up of the
% stationary hazard assessment services in the upper-bounded
% Gutenberg-Richter estimation mode.
%
% The mean return period of magnitude M is the average elapsed time between
% the consecutive earthquakes of magnitude M.
% The mean return periods are calculated for magnitude starting from Md up
% to Mu with step dM.
%
% INPUT:
% t - vector of earthquake occurrence times
% M - vector of earthquake magnitudes
% Md - starting magnitude for return period calculations
% Mu - ending magnitude for return period calculations
% dM - magnitude step for return period calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
% Mmax - upper limit of magnitude distribution
% OUTPUT:
% m - vector of independent variable (magnitude) m=(Md:dM:Mu)
% T - vector od mean return periods of the same length as m
%
% 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, <http://www.gnu.org/licenses/>.
%
function [m,T]=Ret_periodGRT(Md,Mu,dM,Mmin,lamb,eps,b,Mmax)
if Md<Mmin; Md=Mmin;end
if Mu>Mmax; Mu=Mmax;end
m=(Md:dM:Mu)';
beta=b*log(10);
T=1/lamb./(1-Cdfgr(m,beta,Mmin-eps/2,Mmax));
end
function [y]=Cdfgr(t,beta,Mmin,Mmax)
%CDF of the truncated upper-bounded exponential distribution (truncated G-R
% model
% Mmin - catalog completeness level
% Mmax - upper limit of the distribution
% beta - the distribution parameter
% t - vector of magnitudes (independent variable)
% y - CDF vector
mian=(1-exp(-beta*(Mmax-Mmin)));
y=(1-exp(-beta*(t-Mmin)))/mian;
idx=find(y>1);
y(idx)=ones(size(idx));
end

55
src/StationarySeismicHazardAnalysis/Ret_periodGRU.m Normal file
View File

@@ -0,0 +1,55 @@
% [m,T]=Ret_periodGRU(Md,Mu,dM,Mmin,lamb,eps,b)
%
% EVALUATES THE MEAN RETURN PERIOD VALUES USING THE UNLIMITED G-R LED
% MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the unlimited Gutenberg-Richter relation
% leads to the exponential distribution model of magnitude distribution
% from and above the catalog completness level Mmin. The shape parameter of
% this distribution and consequently the G-R b-value are calculated at
% start-up of the stationary hazard assessment services in the
% unlimited Gutenberg-Richter estimation mode.
%
% The mean return period of magnitude M is the average elapsed time between
% the consecutive earthquakes of magnitude M.
% The mean return periods are calculated for magnitude starting from Md up
% to Mu with step dM.
%
%INPUT:
% Md - starting magnitude for return period calculations
% Mu - ending magnitude for return period calculations
% dM - magnitude step for return period calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
%
%OUTPUT:
% m - vector of independent variable (magnitude) m=(Md:dM:Mu)
% T - vector od mean return periods of the same length as m
%
% 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, <http://www.gnu.org/licenses/>.
%
function [m,T]=Ret_periodGRU(Md,Mu,dM,Mmin,lamb,eps,b)
if Md<Mmin; Md=Mmin;end
m=(Md:dM:Mu)';
beta=b*log(10);
T=1/lamb./exp(-beta*(m-Mmin+eps/2));
end

90
src/StationarySeismicHazardAnalysis/Ret_periodNPT.m Normal file
View File

@@ -0,0 +1,90 @@
% [m,T]=Ret_periodNPT(Md,Mu,dM,Mmin,lamb,eps,h,xx,ambd,Mmax)
%
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES THE MEAN
% RETURN PERIOD VALUES FOR THE UPPER-BOUNDED NONPARAMETRIC
% DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution has a hard end point Mmax from the right hand
% side.The estimation makes use of the previously estimated parameters
% namely the mean activity rate lamb, the length of magnitude round-off
% interval, eps, the smoothing factor, h, the background sample, xx, the
% scaling factors for the background sample, ambd, and the end-point of
% magnitude distribution Mmax. The background sample,xx, comprises the
% randomized values of observed magnitude doubled symmetrically with
% respect to the value Mmin-eps/2.
%
% The mean return periods are calculated for magnitude starting from Md up
% to Mu with step dM.
%
% 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
%
% INPUT:
% Md - starting magnitude for return period calculations
% Mu - ending magnitude for return period calculations
% dM - magnitude step for return period calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
% Mmax - upper limit of magnitude distribution
%
% OUTPUT:
% m - vector of independent variable (magnitude) m=(Md:dM:Mu)
% T - vector od mean return periods of the same length as m
%
% 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, <http://www.gnu.org/licenses/>.
%
function [m,T]=Ret_periodNPT(Md,Mu,dM,Mmin,lamb,eps,h,xx,ambd,Mmax)
if Md<Mmin; Md=Mmin;end
if Mu>Mmax; Mu=Mmax;end
m=(Md:dM:Mu)';
n=length(m);
mian=2*(Dystr_npr(Mmax,xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
for i=1:n
CDF_NPT=2*(Dystr_npr(m(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h))/mian;
T(i)=1/lamb./(1-CDF_NPT);
end
T=T';
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

87
src/StationarySeismicHazardAnalysis/Ret_periodNPU.m Normal file
View File

@@ -0,0 +1,87 @@
% [m,T]=Ret_periodNPU(Md,Mu,dM,Mmin,lamb,eps,h,xx,ambd)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES
% THE MEAN RETURN PERIOD VALUES FOR THE UNBOUNDED
% NONPARAMETRIC DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution is unlimited from the right hand side.
% The estimation makes use of the previously estimated parameters of kernel
% estimation, namely the smoothing factor, the background sample and the
% scaling factors for the background sample. The background sample
% - xx comprises the randomized values of observed magnitude doubled
% symmetrically with respect to the value Mmin-eps/2.
%
% The mean return period of magnitude M is the average
% elapsed time between the consecutive earthquakes of magnitude M.
% The mean return periods are calculated for magnitude starting from Md up
% to Mu with step dM.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation fro 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
%
% INPUT:
% Md - starting magnitude for return period calculations
% Mu - ending magnitude for return period calculations
% dM - magnitude step for return period calculations
% Mmin - lower bound of the distribution - catalog completeness level
% lamb - mean activity rate for events M>=Mmin
% eps - length of the round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
%
%OUTPUT:
% m - vector of independent variable (magnitude) m=(Md:dM:Mu)
% T - vector od mean return periods of the same length as m
%
% 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, <http://www.gnu.org/licenses/>.
%
function [m,T]=Ret_periodNPU(Md,Mu,dM,Mmin,lamb,eps,h,xx,ambd)
if Md<Mmin; Md=Mmin;end
m=(Md:dM:Mu)';
n=length(m);
for i=1:n
CDF_NPU=2*(Dystr_npr(m(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
T(i)=1/lamb./(1-CDF_NPU);
end
T=T';
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

241
src/StationarySeismicHazardAnalysis/TruncGR.m Normal file
View File

@@ -0,0 +1,241 @@
%
% [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err]=TruncGR(t,M,iop,Mmin)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% PART OF THE SAMPLE COMPRISING EVENTS >=Mmin (COMPLETE PART),
%THE ROUND-OFF ERROR OF MAGNITUDE, THE GUTENBERG-RICHTER B-VALUE
% AND THE UPPER BOUND OF MAGNITUDE DISTRIBUTION USING THE UPPER-BOUNDED
% G-R LED MAGNITUDE DISTRIBUTION MODEL
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the upper-bounded Gutenberg-Richter
% relation leads to the upper truncated exponential distribution to model
% magnitude distribution from and above the catalog completness level
% Mmin. The shape parameter of this distribution and consequently the G-R
% b-value is estimated by maximum likelihood method (Aki-Utsu procedure).
% 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)).
% The mean activity rate, lamb, is the number of events >=Mmin into the
% length of the period in which they occurred. Upon the value of the input
% parameter, iop, the used unit of time can be either day ot month or year.
% The round-off interval length - eps is the least non-zero difference
% between sample data or 0.1 if the least difference is greater than 0.1.
%
% REFERENCES:
%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 magnitudes from a user selected catalog
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
% iop=2 - 'year'
% Mmin - catalog completeness level. Can take any value from [min(M), max(M)].
%
% 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 15 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 the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
% 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.
%
% 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, <http://www.gnu.org/licenses/>.
%
function [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err]=TruncGR(t,M,iop,Mmin)
n=length(M);
lamb_err=0;
t1=t(1);
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
if iop==0
lamb_all=n/round(t(n)-t(1));
lamb=nn/round(t2-t1);
unit='day';
elseif iop==1
lamb_all=30*n/(t(n)-t(1));
lamb=30*nn/(t2-t1);
unit='month';
else
lamb_all=365*n/(t(n)-t(1));
lamb=365*nn/(t2-t1);
unit='year';
end
if nn<15
eps=0;b=0;Mmax=0;err=0;
lamb_err=1;
return;
end
eps=magn_accur(M);
xx=M(M>=Mmin);
clear x;
nn=length(xx);
Max_obs=max(xx);
beta0=0;
Mmax1=Max_obs;
for i=1:50,
beta=fzero(@bet_est,[0.05,4.0],[],mean(xx),Mmin-eps/2,Mmax1);
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))
err=0;
break;
end
Mmax1=Mmax;
beta0=beta;
end
if i==50;
err=1.0;
Mmax=2*xx(1)-xx(2);
beta=fzero(@bet_est,1.0,[],mean(xx),Mmin-eps/2,Mmax);
end
b=beta/log(10);
clear xx
end
function [zero]=bet_est(b,ms,Mmin,Mmax)
%First derivative of the log likelihood function of the upper-bounded
% exponential distribution (truncated GR model)
% b - parameter of the distribution 'beta'
% ms - mean of the observed magnitudes
% Mmin - catalog completeness level
% Mmax - upper limit of the distribution
M_max_min=Mmax-Mmin;
e_m=exp(-b*M_max_min);
zero=1/b-ms+Mmin-M_max_min*e_m/(1-e_m);
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
function [y]=f_podc(z,n,beta,Mmin,Mmax)
% Integrated function for Mmax estimation. Truncated GR model
% z - column vector of independent variable
% n - the size of 'z'
% beta - the distribution parameter
% Mmin - the catalog completeness level
% Mmax - the upper limit of the distribution
y=Cdfgr(z,beta,Mmin,Mmax).^n;
end
function [y]=Cdfgr(t,beta,Mmin,Mmax)
%CDF of the truncated upper-bounded exponential distribution (truncated G-R
% model
% Mmin - catalog completeness level
% Mmax - upper limit of the distribution
% beta - the distribution parameter
% t - vector of magnitudes (independent variable)
% y - CDF vector
mian=(1-exp(-beta*(Mmax-Mmin)));
y=(1-exp(-beta*(t-Mmin)))/mian;
idx=find(y>1);
y(idx)=ones(size(idx));
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

149
src/StationarySeismicHazardAnalysis/UnlimitGR.m Normal file
View File

@@ -0,0 +1,149 @@
% [lamb_all,lamb,lamb_err,unit,eps,b]=UnlimitGR(t,M,iop,Mmin)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% PART OF THE SAMPLE COMPRISING EVENTS >=Mmin (COMPLETE PART),
%THE ROUND-OFF ERROR OF MAGNITUDE AND THE GUTENBERG-RICHTER B-VALUE
%USING THE UNLIMITED G-R LED MAGNITUDE DISTRIBUTION MODEL
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the unlimited Gutenberg-Richter relation
% leads to the exponential distribution model of magnitude distribution
% from and above the catalog completness level Mmin. The shape parameter of
% this distribution and consequently the G-R b-value is estimated by
% maximum likelihood method (Aki-Utsu procedure).
% The mean activity rate, lamb, is the number of events >=Mmin into the
% length of the period in which they occurred. Upon the value of the input
% parameter, iop, the used unit of time can be either day ot month or year.
% The round-off interval length - eps is either the least non-zero difference
% between sample data or 0.1 if the least difference is greater than 0.1.
%
% INPUT:
% t - vector of earthquake occurrence times
% M - vector of magnitudes from a user selected catalog
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
% iop=2 - 'year'
% Mmin - catalog completeness level. Can take any value from [min(M), max(M)].
%
% 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 7 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 the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
%
% 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 <http://www.gnu.org/licenses/>.
%
function [lamb_all,lamb,lamb_err,unit,eps,b]=UnlimitGR(t,M,iop,Mmin)
if isempty(t) || numel(t)<3 || isempty(M(M>=Mmin))
t=[1 2];M=[1 2]; end
lamb_err=0;
n=length(M);
t1=t(1);
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
[NM,unit]=time_diff(t(1),t(n),iop);
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);
lamb=nn/NM;
if nn<7
eps=0;b=0;
lamb_err=1;
return;
end
eps=magn_accur(M);
xx=M(M>=Mmin);
clear x;
beta=1/(mean(xx)-Mmin+eps/2);
b=beta/log(10);
clear xx
end
function [NM,unit]=time_diff(t1,t2,iop)
% 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 [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

60
src/StationarySeismicHazardAnalysis/dist_GRT.m Normal file
View File

@@ -0,0 +1,60 @@
% [m, PDF_GRT, CDF_GRT]=dist_GRT(Md,Mu,dM,Mmin,eps,b,Mmax)
%
% EVALUATES THE DENSITY AND CUMULATIVE DISTRIBUTION FUNCTIONS OF MAGNITUDE
% UNDER THE UPPER-BOUNDED G-R LED MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the upper-bounded Gutenberg-Richter
% relation leads to the upper truncated exponential distribution to model
% magnitude distribution from and above the catalog completness level
% Mmin. The shape parameter of this distribution, consequently the G-R
% b-value and the end-point of the distribution Mmax are calculated at
% start-up of the stationary hazard assessment services in the
% upper-bounded Gutenberg-Richter estimation mode.
%
% The distribution function values are calculated for magnitude starting
% from Md up to Mu with step dM.
%
%INPUT:
% Md - starting magnitude for distribution functions calculations
% Mu - ending magnitude for distribution functions calculations
% dM - magnitude step for distribution functions calculations
% Mmin - lower bound of the distribution - catalog completeness level
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
% Mmax - upper limit of magnitude distribution
%
%OUTPUT:
% m - vector of the independent variable (magnitude) m=(Md:dM:Mu)
% PDF_GRT - PDF vector of the same length as m
% CDF_GRT - CDF vector of the same length as m
%
% 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, <http://www.gnu.org/licenses/>.
%
function [m, PDF_GRT, CDF_GRT]=dist_GRT(Md,Mu,dM,Mmin,eps,b,Mmax)
m=(Md:dM:Mu)';
beta=b*log(10);
mian=(1-exp(-beta*(Mmax-Mmin+eps/2)));
PDF_GRT=beta*exp(-beta*(m-Mmin+eps/2))/mian;
CDF_GRT=(1-exp(-beta*(m-Mmin+eps/2)))/mian;
idx=find(CDF_GRT<0);
PDF_GRT(idx)=zeros(size(idx));CDF_GRT(idx)=zeros(size(idx));
idx=find(CDF_GRT>1);
PDF_GRT(idx)=zeros(size(idx));CDF_GRT(idx)=ones(size(idx));
end

57
src/StationarySeismicHazardAnalysis/dist_GRU.m Normal file
View File

@@ -0,0 +1,57 @@
% [m, PDF_GRU, CDF_GRU]=dist_GRU(Md,Mu,dM,Mmin,eps,b)
%
% EVALUATES THE DENSITY AND CUMULATIVE DISTRIBUTION FUNCTIONS OF MAGNITUDE
% UNDER THE UNLIMITED G-R LED MAGNITUDE DISTRIBUTION MODEL.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The assumption on the unlimited Gutenberg-Richter relation
% leads to the exponential distribution model of magnitude distribution
% from and above the catalog completness level Mmin. The shape parameter of
% this distribution and consequently the G-R b-value are calculated at
% start-up of the stationary hazard assessment services in the
% unlimited Gutenberg-Richter estimation mode.
%
% The distribution function values are calculated for magnitude starting
% from Md up to Mu with step dM.
%
%INPUT:
% Md - starting magnitude for distribution functions calculations
% Mu - ending magnitude for distribution functions calculations
% dM - magnitude step for distribution functions calculations
% Mmin - lower bound of the distribution - catalog completeness level
% eps - length of the round-off interval of magnitudes.
% b - Gutenberg-Richter b-value
%
%OUTPUT:
% m - vector of the independent variable (magnitude) m=(Md:dM:Mu)
% PDF_GRT - PDF vector of the same length as m
% CDF_GRT - CDF vector of the same length as m
%
%
% 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, <http://www.gnu.org/licenses/>.
%
function [m, PDF_GRU, CDF_GRU]=dist_GRU(Md,Mu,dM,Mmin,eps,b)
m=(Md:dM:Mu)';
beta=b*log(10);
PDF_GRU=beta*exp(-beta*(m-Mmin+eps/2));
CDF_GRU=1-exp(-beta*(m-Mmin+eps/2));
idx=find(CDF_GRU<0);
PDF_GRU(idx)=zeros(size(idx));CDF_GRU(idx)=zeros(size(idx));
idx=find(CDF_GRU>1);
PDF_GRU(idx)=zeros(size(idx));CDF_GRU(idx)=ones(size(idx));
end

110
src/StationarySeismicHazardAnalysis/dist_NPT.m Normal file
View File

@@ -0,0 +1,110 @@
% [m,PDF_NPT,CDF_NPT]=dist_NPT(Md,Mu,dM,Mmin,eps,h,xx,ambd,Mmax)
%
% USING THE NONPARAMETRIC ADAPTATIVE KERNEL ESTIMATORS EVALUATES THE DENSITY
% AND CUMULATIVE DISTRIBUTION FUNCTIONS FOR THE UPPER-BOUNDED MAGNITUDE
% DISTRIBUTION.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution has a hard end point Mmax from the right hand
% side.The estimation makes use of the previously estimated parameters
% namely the mean activity rate lamb, the length of magnitude round-off
% interval, eps, the smoothing factor, h, the background sample, xx, the
% scaling factors for the background sample, ambd, and the end-point of
% magnitude distribution Mmax. The background sample,xx, comprises the
% randomized values of observed magnitude doubled symmetrically with
% respect to the value Mmin-eps/2.
%
% 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
%
%INPUT:
% Md - starting magnitude for distribution functions calculations
% Mu - ending magnitude for distribution functions calculations
% dM - magnitude step for distribution functions calculations
% Mmin - lower bound of the distribution - catalog completeness level
% eps - length of round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
% Mmax - upper limit of magnitude distribution
%
% OUTPUT:
% m - vector of the independent variable (magnitude)
% PDF_NPT - PDF vector
% CDF_NPT - CDF vector
%
% 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 , <http://www.gnu.org/licenses/>.
%
function [m,PDF_NPT,CDF_NPT]=dist_NPT(Md,Mu,dM,Mmin,eps,h,xx,ambd,Mmax)
m=(Md:dM:Mu)';
nn=length(m);
mian=2*(Dystr_npr(Mmax,xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
for i=1:nn
if m(i)<Mmin-eps/2
PDF_NPT(i)=0;CDF_NPT(i)=0;
elseif m(i)>Mmax
PDF_NPT(i)=0;CDF_NPT(i)=1;
else
PDF_NPT(i)=dens_npr1(m(i),xx,ambd,h,Mmin-eps/2)/mian;
CDF_NPT(i)=2*(Dystr_npr(m(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h))/mian;
end
end
PDF_NPT=PDF_NPT';CDF_NPT=CDF_NPT';
end
function [gau]=dens_npr1(y,x,ambd,h,x1)
%Nonparametric adaptive density for a variable from the interval [x1,inf)
% x - the sample data doubled and sorted in the ascending order.
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
c=sqrt(2*pi);
if y<x1
gau=0;
else
gau=2*sum(exp(-0.5*(((y-x)./ambd')./h).^2)./ambd')/c/n/h;
end
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end

109
src/StationarySeismicHazardAnalysis/dist_NPU.m Normal file
View File

@@ -0,0 +1,109 @@
% [m, PDF_NPU, CDF_NPU]=dist_NPU(Md,Mu,dM,Mmin,eps,h,xx,ambd)
%
% USING THE NONPARAMETRIC ADAPTATIVE KERNEL ESTIMATORS EVALUATES THE DENSITY
% AND CUMULATIVE DISTRIBUTION FUNCTIONS FOR THE UNLIMITED MAGNITUDE
% DISTRIBUTION.
%
% AUTHOR: Stanislaw. Lasocki, Institute of Geophysics Polish Academy of
% Sciences, Warsaw, Poland
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. It is assumed that
% the magnitude distribution is unlimited from the right hand side.
% The estimation makes use of the previously estimated parameters of kernel
% estimation, namely the smoothing factor, the background sample and the
% scaling factors for the background sample. The background sample
% - xx comprises the randomized values of observed magnitude doubled
% symmetrically with respect to the value Mmin-eps/2
%
% The distribution function values are calculated for magnitude starting
% from Md up to Mu with step dM.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation fro 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
%
%INPUT:
% Md - starting magnitude for distribution functions calculations
% Mu - ending magnitude for distribution functions calculations
% dM - magnitude step for distribution functions calculations
% Mmin - lower bound of the distribution - catalog completeness level
% eps - length of round-off interval of magnitudes.
% h - kernel smoothing factor.
% xx - the background sample
% ambd - the weigthing factors for the adaptive kernel
%
%
%OUTPUT
% m - vector of the independent variable (magnitude) m=(Md:dM:Mu)
% PDF_NPU - PDF vector of the same length as m
% CDF_NPU - CDF vector of the same length as m
%
% 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 , <http://www.gnu.org/licenses/>.
%
function [m, PDF_NPU, CDF_NPU]=dist_NPU(Md,Mu,dM,Mmin,eps,h,xx,ambd)
m=(Md:dM:Mu)';
nn=length(m);
for i=1:nn
if m(i)>=Mmin-eps/2
PDF_NPU(i)=dens_npr1(m(i),xx,ambd,h,Mmin-eps/2);
CDF_NPU(i)=2*(Dystr_npr(m(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
else
PDF_NPU(i)=0;
CDF_NPU(i)=0;
end
end
PDF_NPU=PDF_NPU';CDF_NPU=CDF_NPU';
end
function [gau]=dens_npr1(y,x,ambd,h,x1)
%Nonparametric adaptive density for a variable from the interval [x1,inf)
% x - the sample data doubled and sorted in the ascending order.
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
c=sqrt(2*pi);
if y<x1
gau=0;
else
gau=2*sum(exp(-0.5*(((y-x)./ambd')./h).^2)./ambd')/c/n/h;
end
end
function [Fgau]=Dystr_npr(y,x,ambd,h)
%Nonparametric adaptive cumulative distribution for a variable from the
%interval (-inf,inf)
% x - the sample data
% ambd - the local scaling factors for the adaptive estimation
% h - the optimal smoothing factor
% y - the value of random variable X for which the density is calculated
% gau - the density value f(y)
n=length(x);
Fgau=sum(normcdf(((y-x)./ambd')./h))/n;
end