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

SERA Toolbox1 and Toolbox2 standalone versions

Browse Source
This commit is contained in:
Konstantinos Leptokaropoulos
2019-07-05 10:31:31 +02:00
parent 6f5c52a565
commit 758751a7b0
71 changed files with 54733 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

701
Magnitude_Complexity_TOOLBOX_D23_2/Multimodality_Test_VERSION_1/MM_MB_V1_8.m Normal file
View File

@@ -0,0 +1,701 @@
% FUNCTION: MM_MB
% VERSION: [Interactive Hybrid Standalone Version] V1.8
% COMPATIBLE with Matlab version 2017b or later
% TOOLBOX: "Magnitude Complexity Toolbox" within SERA Project
% DOCUMENT: "READ_ME_App_2B_v1_Description_MM_MB.docx"
%--------------------------------------------------------------------------------------------------------------
%% EXAMPLE TO RUN:
% x=exprnd(log10(exp(1)),1000,1);
% [n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps,gau,gau_b,poch,poch2,zer1,zer2,x1,x2]...
% =MM_MB_V1_8(x,0.2,500,500,0.01,0.001,'Efron','Silverman') % for Interactivity OFF
% [n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps,gau,gau_b,poch,poch2,zer1,zer2,x1,x2]...
% =MM_MB_V1_8(x); % for Interactivity ON
%% ----------------------------------------------------------------------------------------------------------
% Test performed for Magnitude Distribution of a given dataset checking
% whether PDF demonstrates multi-modes/ multi-bumps
% --------------------------------------------------------------------------------------------------------
% INPUT: THE CURRENT VERSION USES AS INPUT ANY MAGNITUDE VECTOR
% (Appropriate for standalone use - function mode)
% --------------------------------------------------------------------------------------------------------
% OVERVIEW: This Application is a Matlab function which performs testing
% of hypotheses of 1) multimodality and 2) existence of multi-bumps in a given
% magnitude distribution.
% --------------------------------------------------------------------------------------------------------
% AUTHORS: K. Leptokaropoulos and P. Urban,
% last updated: 03/2019, within SERA PROJECT, EU Horizon 2020 R&I
% programme under grant agreement No.730900
% CURRENT VERSION: v1.8 **** [INTERACTIVE HYBRID STANDALONE VERSION!!]
% ----- THIS IS a dual-mode version: If all input arguments are set, then the
% ----- Application operates as a function. However, if only the input vector
% ----- is introduced, the application switch to interactive mode.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PLEASE refer to the accompanying document:
% "READ_ME_App_2B_v1_Description_MM_MB.docx"
% for a description of the Application and its requirements.
% --------------------------------------------------------------------------------------------------------
% DESCRIPTION:
% The function studies the magniutude distribution complexity by means of the
% Mulitmodality Test (Silverman, 1986; Efron and Tibshirani, 1993). Two null
% hypotheses (H0s) are tested:
% H01 - multimodality: The magnitude PDF is unimodal
% H02 - multi-bump: The magnitude PDF has one bump to the right of the
% mode.
% A mode is a local maximum of probability density and a bump is an interval
% [a,b] such that the probability density is concave over [a,b] but not over any
% larger interval (Silverman, 1986). The importance of modes and bumps relies
% on the fact that multiple occurrences of these features in a PDF indicate, for
% most standard densities, a mixing of components (e.g. Cox, 1966).
% --------------------------------------------------------------------------------------------------------
% INPUT: If the User Provides only the magnitude vector as input, then the
% function enables a GUI allowing the User to interactivelly set the rest
% of the input parameters by the User. Otherwise the User has to define
% the 7 input parameters in addition to the input magnitude vector.
% PLEASE REFER TO "MM_MB.doc" FOR A COMPREHENSIVE
% PARAMETER DESCRIPTION
% Input Parameters Overview:
% --- M: A Time-Series (e.g. Magnitude) Vector. The program
% takes as input any matlab vector. The input data can be
% uploaded by the User from an ASCII file (e.g. the file
% "test_vector.txt", located in the "Sample_Data"). Such file
% should contain a vector (raw or column) of the Data the
% User wishes to process. The User is afterwards requested
% to enter parameters values:
% --- Mc: Corresponds to the catalog completeness threshold.
% --- m: number of points to divide magnitude sample, default 100,
% recommended 100-1000.
% % INPUTS for MODE and BUMP HUNTING
% - n_boot: Number of bootstrap iterations for both MM and MB
% (default: 100)
% - delta_h: Smoothing parameter step for successive trials in defining
% the critical h for MM and MB process (h-critical accuracy)
% (default: 0.001)
% - h : initial value of the smoothing factor to apply in defining the
% critical h for the MM testing process (default:0.01)
% - MMmeth: Method for multimodality Testing. Possible values:
% 'Efron', 'Silverman'
% - MBmeth: Method for multibumps Testing. Possible values:
% 'Silverman', 'Efron'
% NOTE: The initial value for MB process is taken equal to the
% critical value defined in the MM process minus its
% accuracy (i.e. hcrit_Modes-delta_h, see OUTPUT below)
% ---------------------------------------------------------------------------------------------------------
% OUTPUT:
% - Output Report with data and parameters used
% ('Output_MM_MB.txt file')
% - Output Parameters:
% % DOUBLES:
% * n - Number of observations used
% * bval- b-value of the G-R law
% * Rmodes - The estimated significance of null hypothesis (H01)
% * hcrit_modes - estimated critical smoothness parameter for MM test
% * Rbumps - The estimated significance of null hypothesis (H02)
% * hcrit_bumps - estimated critical smoothness parameter for MB test
% * zer1 - point where 1st derivative is zero (extremum)
% * zer2 - point(s) where 2nd derivative is zero (inflection points)
% [NOTE: zer2 can be a vector]
% % VECTORS:
% * gau: Data vector PDF estimated for h critical from Multimodality
% * gau_b: Data vector PDF estimated for h critical from Multibumps
% * poch: 1st derivative of Data vector PDF for h critical from Multibumps
% * poch2: 2nd derivative of Data vector PDF for h critical from Multibumps
% * x1,x2: minimum and maximum Data points after randomization
% (parameters needed for plotting)
% - Output Figures:
%
% ---------------------------------------------------------------------------------------------------------
% REFERENCES:
% -- Cox, D. R., (1966), Notes on the analysis of mixed frequency distributions.
% Br. J. Math. Stat. Psychol. 19, 39-47, doi.org/10.1111/j.2044-8317.1966.tb00353.x.
% -- Efron, B., and Tibshirani R.J. (1993) An Introduction to the Bootstrap,
% CRC Press, Boca Raton, Fla.
% -- Lasocki S. and E. E, Papadimitriou (2006), "Magnitude distribution
% complexity revealed in seismicity from Greece", J. Geophys. Res.,
% 111, B11309, doi:10.1029/2005JB003794.
% -- Silverman, B. W., (1986), Density estimation for statistics and data analysis,
% CRC press, Boa Raton, Fla.
% ---------------------------------------------------------------------------------------------------------
% LICENSE
% 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.
% ---------------------------------------------------------------------------------------------------------
%% EXAMPLE TO RUN:
% x=exprnd(log10(exp(1)),1000,1);
% [n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps,gau,gau_b,poch,poch2,zer1,zer2,x1,x2]...
% =MM_MB_V1_8(x,0.2,500,500,0.01,0.001,'Efron','Silverman') % for Interactivity OFF
% [n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps,gau,gau_b,poch,poch2,zer1,zer2,x1,x2]...
% =MM_MB_V1_8(x); % for Interactivity ON
%% ----------------------------------------------------------------------------------------------------------
function [n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps,gau,gau_b,poch,poch2,zer1,zer2,x1,x2]...
=MM_MB_V1_8(M,Mc,m,n_boot,h,delta_h,MMmeth,MBmeth)
mkdir Outputs_MM_MB
% define round-off interval and select data above Mc
sm=sort(M);dm=sm(2:length(sm))-sm(1:length(sm)-1);
dm=dm(dm>0);eps_M=min(dm);
%% --------------------------------------- INTERACTIVE MODE ----------------------------------------
if nargin==1
[Mi,Mmin]=FiltMcVector(M,eps_M);
M=Mi(Mi>=Mmin-eps_M/2);
m=dialog1('number of points to divide sample',{'100'});
Rp=round(-log10(eps_M)); %% check the RP parameter
[beta]=beta_AK_UT_Mbin (Mmin,M,Rp) ;
% Randomize magnitude
m_corr(:,1) = korekta(M,Mmin,eps_M,beta);
% SETTING PARAMETERS
n=length(m_corr);
x1=floor(min(m_corr)*10)/10; %nieco mniej ni<EFBFBD> min(m) zaokr<EFBFBD>glone do jednej dziesi<EFBFBD>tej
x2=ceil(max(m_corr)*10)/10; %nieco wi<EFBFBD>cej ni<EFBFBD> max(m) zaokr<EFBFBD>glone do jednej dziesi<EFBFBD>tej
% MULTIMODALITY TESTING
% Give number of bootstrap attempts and h values
[n_boot,delta_h,h]=dialog2('Multimodality',{'100','0.001','0.01'});
n_boot1=n_boot;
%delta_h=0.1;%input('czynnik wygladzajacy '); %0.1
%h=0.01; %0.001
[hcrit,gau]=critical_smoothing(m_corr,n,delta_h,x1,x2,m,h);
hcrit;
% Select Method of Multimodality testing 'Silverman' or 'Efron'
clc;
disp('Please Select Method for Multimodality Testing: ');
disp('(1 - for Efron, 2 - for Silverman)');
MM=input(' ');
% ---------testing plot-----------
plot(1:m,gau);title('Data density for critical smoothing factor (h_c_r_i_t)');
disp('press any key ...')
pause
close all
% ---------testing plot-----------
% Run Multimodality testing
switch MM
case 1
[R] = test_multimodality_e (m_corr,n,x1,x2,m,hcrit,n_boot1);
case 2
[R] = test_multimodality (m_corr,n,x1,x2,m,hcrit,n_boot1);
end
% ------------------------------------------------------------------------
% [zer1]=zera_1st(m_corr,n,x1,x2,m,hcrit-delta_h);
% ------------------------------------------------------------------------
% BUMP HUNTING
% Give number of bootstrap attempts and h values
[n_boot,delta_h,h]=dialog2('Multibumps',{'100','0.001',num2str(hcrit)});
n_boot2=n_boot;
%delta_h=0.1;
[hcrit_bumps,gau_b,poch,poch2]=critical_smoothing_bumps(m_corr,n,delta_h,x1,x2,m,hcrit-delta_h);
hcrit_bumps;
% Select Method for testing the bump size 'Silverman' or 'Efron'
disp(' ');disp(' ');disp(' ');
disp('Please Select Method for Bump Testing: ');
disp('(1 - for Silverman, 2 - for Efron)');
BM=input(' ');
% Run Bump hunting
switch BM
case 1
[Rb] = bump_hunt (m_corr,n,x1,x2,m,hcrit_bumps,n_boot2);
case 2
[Rb] = bump_hunt_e (m_corr,n,x1,x2,m,hcrit_bumps,n_boot2);
end
%% -------------------------------------- FUNCTION MODE ----------------------------------------
elseif nargin>1
% INPUT PARAMETERS ARE SPECIFIED AS FUNCTION ARGUMENTS
MM=M;Mmin=Mc;M=M(M>=Mmin-eps_M/2);
m;n_boot;h;delta_h;MMmeth;MBmeth;
Rp=round(-log10(eps_M)); %% check the RP parameter
[beta]=beta_AK_UT_Mbin (Mmin,M,Rp); % beta value
% Randomize magnitude
m_corr(:,1) = korekta(M,Mmin,eps_M,beta);
% SETTING PARAMETERS
n=length(m_corr);
x1=floor(min(m_corr)*10)/10; %nieco mniej ni<EFBFBD> min(m) zaokr<EFBFBD>glone do jednej dziesi<EFBFBD>tej
x2=ceil(max(m_corr)*10)/10; %nieco wi<EFBFBD>cej ni<EFBFBD> max(m) zaokr<EFBFBD>glone do jednej dziesi<EFBFBD>tej
% MULTIMODALITY TESTING
[hcrit,gau]=critical_smoothing(m_corr,n,delta_h,x1,x2,m,h);
switch MMmeth
case 'Efron'
[R] = test_multimodality_e (m_corr,n,x1,x2,m,hcrit,n_boot);
case 'Silverman'
[R] = test_multimodality (m_corr,n,x1,x2,m,hcrit,n_boot);
end
% BUMP HUNTING
[hcrit_bumps,gau_b,poch,poch2]=critical_smoothing_bumps(m_corr,n,delta_h,x1,x2,m,hcrit-delta_h);
switch MBmeth
case 'Silverman'
[Rb] = bump_hunt (m_corr,n,x1,x2,m,hcrit_bumps,n_boot);
case 'Efron'
[Rb] = bump_hunt_e (m_corr,n,x1,x2,m,hcrit_bumps,n_boot);
end
end
bval=beta/log(10);
% Display the Results:
re=table(n, bval, R,hcrit,Rb,hcrit_bumps)
disp(['Significance of Ho(1) that input data distribution is unimodal is: ', num2str(R)])
disp(['Significance of Ho(2) that input data distribution has no more '])
disp(['than one bumps to the right of the global maximum is: ', num2str(Rb)])
[zer1]=zera_1st(m_corr,n,x1,x2,m,hcrit_bumps-delta_h); % use the same h in both cases???
[zer2]=zera_2nd(m_corr,n,x1,x2,m,hcrit_bumps-delta_h);
%% -------------------------- PLOTTING AND SAVING RESULTS ----------------------------
% ----- plotting magnitude density and its derivatives ------
xp=linspace(x1,x2,m);
%st=range([x1 x2])/(m-1);xp=x1:st:x2;
figure('rend','painters','pos',[800 100 600 900]);hold on
subplot(3,1,1);
plot(xp,gau_b,'LineWidth',2);hold on;
for i=1:length(zer1)
Le1=plot([zer1(i) zer1(i)],[0 10],'k--','LineWidth',1);
end
for i=1:length(zer2)
Le2=plot([zer2(i) zer2(i)],[0 10],'r--','LineWidth',1);
end
xlim([min(M)-0.1 max(M)+0.1]);ylim([0 max(gau_b)+0.01])
title(['Data density for h_c_r_i_t ';' (critical smoothing factor) '],'FontSize',14);%axis square
legend([Le1,Le2],{'extremum','Inflection Point(s)'});ylabel('PDF','FontSize',14)
subplot(3,1,2);
plot(xp,poch,'LineWidth',2);ylabel('1^s^t derivative','FontSize',14);hold on;%axis square;
plot([xp(1) xp(length(xp))],[0 0],'--','LineWidth',1)
for i=1:length(zer1)
plot([zer1(i) zer1(i)],[min(poch)-0.01 max(poch)+0.01],'k--','LineWidth',1)
end
for i=1:length(zer2)
plot([zer2(i) zer2(i)],[min(poch)-0.01 max(poch)+0.01],'r--','LineWidth',1)
end
xlim([min(M)-0.1 max(M)+0.1]);ylim([min(poch)-0.01 max(poch)+0.01])
subplot(3,1,3);
plot(xp,poch2,'LineWidth',2);ylabel('2^n^d derivative','FontSize',14);%axis square
hold on; plot([xp(1) xp(length(xp))],[0 0],'--','LineWidth',1)
for i=1:length(zer1)
plot([zer1(i) zer1(i)],[min(poch2)-0.01 max(poch2)+0.01],'k--','LineWidth',1)
end
for i=1:length(zer2)
plot([zer2(i) zer2(i)],[min(poch2)-0.01 max(poch2)+0.01],'r--','LineWidth',1)
end
xlim([min(M)-0.1 max(M)+0.1]);ylim([min(poch2)-0.01 max(poch2)+0.01])
xlabel('Data','FontSize',14);
% Save Outputs
cd Outputs_MM_MB
Rmodes=R;Rbumps=Rb;hcrit_modes=hcrit;
SaveOuts(eps_M,Mmin,m,n_boot,n_boot,n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps)
saveas(gcf,'Multimodality_Output.jpg')
cd ../
%% ******************************************************************
% *************************** FUNCTIONS ***************************
% ****-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-****
%% Select Mc and filter parameters for M>=Mc
function [Cmag,Mc]=FiltMcVector(magnitude_vector,EPS1)
Cmag=magnitude_vector;
%sm=sort(Cmag); sm1=sm(2:length(sm))-sm(1:length(sm)-1); EPS1=min(sm1(sm1>0));
%EPS=dialog1('magnitude round-off interval',{num2str(EPS1)});
ar=min(Cmag):0.1:max(Cmag);
fig_Mc=figure;histogram(Cmag,length(ar));set(gca,'YScale','log')
title('Please Select Data Cutoff','FontSize',14);xlabel('M','FontSize',14),ylabel('Log_1_0N','FontSize',14)
% Select events above Mc
[Mc,N]=ginput(1);Mc=floor(Mc/EPS1)*EPS1,close(fig_Mc);
%Cmag=Cmag(Cmag>=Mc);
end
%% --------------------------------------------------------------------------------------
%% --------------------------------------------------------------------------------------
% finds a field defined by a certain (string) name
function [id] = findfield( catalog,field )
id=0;
j=1;
while j <= size(catalog,2) && id==0
if (strcmp(catalog(j).field,field)==1)
id=j;
end
j=j+1;
end
end
%% --------------------------------------------------------------------------------------
function [ou]=dialog1(name,defaultanswer)
prompt=['\fontsize{12} Please enter ',name, ':'];
prompt={prompt};
numlines=1; opts.Interpreter='tex';
ou=inputdlg(prompt,'Parameter Setting',numlines,defaultanswer,opts);ou=str2num(ou{1});
end
%% --------------------------------------------------------------------------------------
function [ou1, ou2, ou3]=dialog2(titl,defaultanswer)
prompt={'\fontsize{12} Please enter number of bootstrap iterations:','\fontsize{12} Please enter \Deltah (h_c_r_i_t accuracy):','\fontsize{12} Please enter h (smoothing factor):'};
name=['Parameters for ',titl, ':'];
numlines=1;
opts.Interpreter='tex';
answer=inputdlg(prompt,name,numlines,defaultanswer,opts);
ou1=answer(1);ou2=answer(2);ou3=answer(3);
ou1=str2num(ou1{1});ou2=str2num(ou2{1});ou3=str2num(ou3{1});
end
%% --------------------------------------------------------------------------------------
%
function [beta]=beta_AK_UT_Mbin(Mmin,m,Rp)
%
% m - magnitude vector
% Mmin - completeness magitude threshold
% beta - beta value. b(G-R)=beta/log(10)
% Rp - Rounding precision, (1 - one decimal, 2 - two decimals, etc)
beta=1/(mean(m)-(Mmin-0.5*10^(-Rp)));
end
%% --------------------------------------------------------------------------------------
% Magnitude randomization
%
function [m_corr]=korekta(m,Mmin,eps,beta)
%
% m - magnitude vector
% Mmin - completeness magitude threshold
% beta - beta value. b(G-R)=beta/log(10)
% EPS - magnitude round-off interval
%
% m_corr - randomized magnitude vector
%
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
%% --------------------------------------------------------------------------------------
% Finding critical smoothing parameter value for which the non parametric
% PDF demonstrates one mode in the selected interval
%
function [hcrit,gau]=critical_smoothing(mm,n,delta_h,x1,x2,m,h)
%
x=sort(mm);
p=2.0;
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,m);
%h=input('Initial smoothing factor ');
while p>1
for i=1:m,
poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./h).^2));
end
p=0;
for i=2:m,
if poch(i-1)*poch(i)<0
p=p+1;
end
end
h=h+delta_h;
end
hcrit=h-delta_h;
for i=1:m,
gau(i)=sum(exp(-0.5*((y(i)-x)./hcrit).^2))/c/n/hcrit;
end
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
%
function [R]=test_multimodality(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_wielom=0;
for j=1:mm,
ind=unidrnd(n,n,1);
eps=normrnd(0,1,n,1);
z=x(ind)+h_crit.*eps;
z=sort(z);
poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
p=0;
for i=2:n_point,
poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
if poch(i-1)*poch(i)<0
p=p+1;
end
end
if p>1
no_wielom=no_wielom+1;
end
end
R=no_wielom/mm;
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
% EFRON
%
function [R]=test_multimodality_e(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
s2=sqrt(1+h_crit^2/var(x));
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_wielom=0;
for j=1:mm,
ind=unidrnd(n,n,1);
eps=normrnd(0,1,n,1);
me=mean(x(ind));
z=me+(x(ind)-me+h_crit.*eps)./s2;
z=sort(z);
poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
p=0;
for i=2:n_point,
poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
if poch(i-1)*poch(i)<0
p=p+1;
end
end
if p>1
no_wielom=no_wielom+1;
end
end
R=no_wielom/mm;
end
%% --------------------------------------------------------------------------------------
% Finding critical smoothing parameter value for which the non parametric
% PDF demonstrates one bump in the selected interval
%
function [hcrit,gau,poch,poch2]=critical_smoothing_bumps(mm,n,delta_h,x1,x2,m,h)
%
x=sort(mm);
p=2.0;
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,m);
%h=input('Initial smoothing factor ');
while p>1
for i=1:m,
% poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./h).^2));
poch2(i)=sum((((y(i)-x)./h).^2-1).*exp(-0.5*((y(i)-x)./h).^2));
end
p=0;
for i=2:m,
if poch2(i-1)*poch2(i)<0
p=p+1;
end
end
h=h+delta_h;
end
hcrit=h-delta_h;
for i=1:m,
gau(i)=sum(exp(-0.5*((y(i)-x)/hcrit).^2))/c/n/hcrit;
poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./hcrit).^2));
poch2(i)=sum((((y(i)-x)./hcrit).^2-1).*exp(-0.5*((y(i)-x)./hcrit).^2));
end
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
%
function [R]=bump_hunt(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_bump=0;
for j=1:mm,
ind=unidrnd(n,n,1);
eps=normrnd(0,1,n,1);
z=x(ind)+h_crit.*eps;
z=sort(z);
% poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
poch2(1)=sum((((y(1)-z)./h_crit).^2-1).*exp(-0.5*((y(1)-z)./h_crit).^2));
p=0;
for i=2:n_point,
% poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
poch2(i)=sum((((y(i)-z)./h_crit).^2-1).*exp(-0.5*((y(i)-z)./h_crit).^2));
if poch2(i-1)*poch2(i)<0
p=p+1;
end
end
if p>1
no_bump=no_bump+1;
end
end
R=no_bump/mm;
end
%% --------------------------------------------------------------------------------------
% "SMOOTHED BOOTSTRAP"
% EFRON
%
function [R]=bump_hunt_e(m,n,x1,x2,n_point,h_crit,mm)
%
x=sort(m);
s2=sqrt(1+h_crit^2/var(x));
c=sqrt(2*pi);
%x1=input('Lower limit ');
%x2=input('Upper limit ');
%m=input('No. of points ');
y=linspace(x1,x2,n_point);
%h_crit=input('Critical smoothing factor ');
%mm=input('No of trials ');
no_bump=0;
for j=1:mm,
ind=unidrnd(n,n,1);
eps=normrnd(0,1,n,1);
me=mean(x(ind));
z=me+(x(ind)-me+h_crit.*eps)./s2;
z=sort(z);
% poch(1)=-sum((y(1)-z).*exp(-0.5*((y(1)-z)./h_crit).^2));
poch2(1)=sum((((y(1)-z)./h_crit).^2-1).*exp(-0.5*((y(1)-z)./h_crit).^2));
p=0;
for i=2:n_point,
% poch(i)=-sum((y(i)-z).*exp(-0.5*((y(i)-z)./h_crit).^2));
poch2(i)=sum((((y(i)-z)./h_crit).^2-1).*exp(-0.5*((y(i)-z)./h_crit).^2));
if poch2(i-1)*poch2(i)<0
p=p+1;
end
end
if p>1
no_bump=no_bump+1;
end
end
R=no_bump/mm;
end
%% --------------------------------------------------------------------------------------
% Determination of PDF first derivative zeros estimates in the selected interval
%
function [zer1]=zera_1st(mm,n,x1,x2,m,h)
%
x=sort(mm);
c=sqrt(2*pi);
y=linspace(x1,x2,m);
for i=1:m,
poch(i)=-sum((y(i)-x).*exp(-0.5*((y(i)-x)./h).^2));
end
p=0;
for i=2:m,
if poch(i-1)*poch(i)<0
p=p+1;
zer1(p)=(y(i)+y(i-1))/2; %K 23JAN2019
end
end
end
%% --------------------------------------------------------------------------------------
% Determination of PDF second derivative zeros in the selected interval
%
function [zer2]=zera_2nd(mm,n,x1,x2,m,h)
%
x=sort(mm);
c=sqrt(2*pi);
y=linspace(x1,x2,m);
for i=1:m,
poch2(i)=sum((((y(i)-x)./h).^2-1).*exp(-0.5*((y(i)-x)./h).^2));
end
p1=0;
for i=2:m,
if poch2(i-1)*poch2(i)<0
p1=p1+1;
zer2(p1)=(y(i)+y(i-1))/2; %K 23JAN2019
end
end
end
%% --------------------------------------------------------------------------------------------------------
% --------------------------------------- SAVE OUTPUTS in the report file ---------------------------------------
% Save Outputs
function SaveOuts(EPS,Mmin,m,n_boot1,n_boot2,n,bval,Rmodes,hcrit_modes,Rbumps,hcrit_bumps)
% ---- Save *.txt file with Parameters Report ----
%cd Outputs/
fid=fopen('REPORT_Multimodality.txt','w');
fprintf(fid,['RESULTS from MULTIMODALITY/MULTIBUMP TESTING (created on ', datestr(now),')\n']);
fprintf(fid,'------------------------------------------------------------------------------\n');
fprintf(fid,['<Round-off interval >: ', num2str(EPS),'\n']);
fprintf(fid,['<Number of points to divide the sample >: ', num2str(m),'\n']);
fprintf(fid,['<Completeness Threshold >: ', num2str(Mmin),'\n']);
fprintf(fid,['<Number of events used >: ', num2str(n),'\n']);
fprintf(fid,['<Gutenberg-Richter b-value >: ', num2str(bval,'%5.3f'),'\n']);
fprintf(fid,['<Number of bootstrap iterations (multimodality) >: ', num2str(n_boot1),'\n']);
fprintf(fid,['<Number of bootstrap iterations (multibumps) >: ', num2str(n_boot2),'\n']);
fprintf(fid,'------------------------------------------------------------------------------\n');
fprintf(fid,['<Critical Smoothing Parameter,h (for modes) >: ', num2str(hcrit_modes,'%6.4f'),'\n']);
fprintf(fid,['<Critical Smoothing Parameter,h (for bumps) >: ', num2str(hcrit_bumps,'%6.4f'),'\n']);
fprintf(fid,['<Significance of H01 that input data PDF is unimodal >: ', num2str(Rmodes,'%4.3f'),'\n']);
fprintf(fid,['<Significance of H02 that input data PDF has no more \n']);
fprintf(fid,[' than one bump to the right of the global maximum >: ', num2str(Rbumps,'%4.3f'),'\n']);
fclose(fid);
end
end

BIN
Magnitude_Complexity_TOOLBOX_D23_2/Multimodality_Test_VERSION_1/READ_ME_MM_MB_v1.docx Normal file
View File

Binary file not shown.

880
Magnitude_Complexity_TOOLBOX_D23_2/Multimodality_Test_VERSION_1/Sample_Data/test_vector.txt Normal file
View File

@@ -0,0 +1,880 @@
1.7
0.8
1.8
1.5
1.7
0.8
1.7
2.1
0.9
1.3
1.1
1.9
2.3
1.7
1.9
1.4
1.5
1.6
1.5
2.7
2.9
1.3
1.2
0.9
2.7
1.4
1.2
1.4
1.7
1.5
1.4
1.5
2.4
2.7
1.7
2.1
1.0
1.7
2.3
1.3
1.7
1.2
1.8
1.5
1.7
1.1
1.6
1.3
1.1
1.3
1.2
1.1
2.1
0.9
1.5
1.5
0.9
2.1
1.5
2.0
3.4
2.7
2.1
1.8
2.1
1.4
2.6
0.8
1.5
4.1
1.7
2.0
2.1
1.7
2.1
2.0
1.7
2.2
1.7
2.1
2.2
1.5
3.6
2.0
1.5
1.5
2.0
1.8
1.3
1.5
2.8
1.3
1.5
1.9
1.5
2.4
1.6
0.5
1.8
1.4
1.6
1.5
2.2
1.4
1.5
1.9
2.5
3.2
2.9
1.0
1.4
1.8
2.0
2.2
1.4
1.6
1.8
2.0
1.3
2.8
1.4
1.3
1.1
1.8
1.6
1.1
0.8
1.4
1.0
1.7
1.3
1.5
1.7
3.2
1.9
1.6
2.1
1.6
1.7
1.9
1.6
2.6
1.2
2.3
2.1
2.1
1.8
1.4
1.1
1.9
3.3
1.4
1.6
1.8
1.7
2.4
1.6
1.7
2.2
2.9
2.7
1.3
2.2
1.4
1.9
1.6
1.4
2.0
1.5
1.5
2.1
1.8
3.3
1.5
1.3
1.9
1.3
1.9
3.8
1.7
1.2
2.2
1.7
1.6
2.3
1.6
1.8
2.7
1.5
1.4
1.5
1.6
1.3
1.6
1.1
2.0
1.8
0.8
2.5
1.7
1.9
1.8
3.2
1.1
1.9
2.9
1.1
1.7
1.8
1.6
1.6
1.9
1.4
1.6
1.5
1.7
1.6
1.8
1.3
1.4
0.6
1.4
1.2
1.8
1.7
1.6
1.3
1.6
1.5
2.4
2.0
2.1
2.5
1.8
1.4
2.0
1.1
1.4
2.5
1.5
1.9
1.9
1.6
1.2
1.3
2.8
2.8
2.7
2.4
2.6
2.3
1.0
1.6
1.3
2.0
0.8
1.7
0.7
1.1
1.2
0.6
1.1
0.9
3.1
0.9
1.0
2.0
1.6
1.1
1.0
1.2
2.3
1.5
2.2
1.2
1.6
2.6
1.4
1.3
1.9
1.6
2.2
1.7
2.0
2.4
1.3
1.6
1.8
1.7
1.7
2.1
2.2
2.3
1.8
2.3
1.7
1.4
1.6
2.5
1.3
1.1
1.4
3.0
1.2
1.7
1.7
1.8
2.2
1.7
2.1
2.9
1.8
1.8
2.1
1.7
1.2
2.3
1.2
1.5
1.7
1.8
1.4
1.5
2.7
2.4
1.6
1.9
2.2
1.6
1.6
1.9
1.7
1.8
1.8
2.0
1.0
1.2
1.3
1.6
2.9
1.5
1.3
1.4
1.3
1.7
1.8
1.9
1.9
3.7
1.5
2.0
1.6
1.6
1.5
2.5
4.2
1.6
3.6
1.9
1.8
2.0
1.8
3.0
2.4
1.2
1.5
2.8
2.8
1.7
1.8
2.3
1.5
1.5
1.9
1.9
1.8
1.2
1.2
1.3
2.1
2.0
1.8
1.7
1.6
1.9
1.9
2.0
1.7
1.8
1.2
2.1
0.8
2.2
1.9
1.6
1.0
2.1
2.3
1.6
1.2
1.9
1.7
2.3
1.8
3.3
1.7
2.5
2.0
1.2
1.5
2.5
1.8
2.7
1.2
3.4
1.6
2.4
1.6
2.2
0.6
2.0
1.9
1.6
2.4
1.4
1.3
1.1
2.3
0.5
0.7
0.8
1.8
1.5
1.0
2.3
1.7
0.5
1.8
2.7
2.5
1.5
2.1
5.8
1.5
1.1
1.5
2.4
2.2
1.2
1.9
1.0
2.0
1.2
1.1
1.5
1.9
0.5
2.5
1.6
1.4
1.9
2.5
1.3
2.1
1.6
1.6
1.3
1.7
1.5
2.1
1.6
1.5
3.2
1.2
2.6
1.4
1.3
1.6
1.7
1.4
1.6
1.8
1.5
1.9
0.9
2.6
1.6
1.8
2.1
1.6
1.2
0.8
1.6
1.2
0.7
1.1
3.1
2.4
2.1
2.2
3.0
1.6
1.8
1.5
3.2
1.1
1.4
1.9
1.2
1.9
1.4
2.4
1.8
1.3
1.8
2.3
1.9
1.9
1.5
1.2
1.6
1.5
2.4
1.9
1.5
1.8
1.7
1.8
2.2
1.5
1.6
2.3
1.8
2.7
1.7
2.0
3.0
1.8
2.1
1.5
1.0
1.9
1.7
2.6
2.7
2.0
1.5
1.9
1.7
2.1
1.7
1.3
1.6
2.9
3.1
1.7
2.4
1.3
2.0
2.0
1.7
4.6
2.6
1.5
2.0
1.3
1.4
1.8
1.4
1.3
1.9
1.4
1.7
1.5
1.7
0.9
1.9
1.3
1.2
1.4
1.4
1.3
1.1
1.3
2.1
1.5
1.9
1.9
4.6
1.1
1.0
1.7
1.2
1.7
0.6
1.3
1.7
1.8
1.9
1.6
1.6
1.7
1.8
2.2
1.7
1.4
1.1
1.1
1.5
1.7
1.8
1.3
1.1
1.7
1.6
0.9
1.8
1.8
1.2
1.3
1.6
0.8
1.4
2.2
1.8
1.5
1.9
2.0
1.7
1.6
1.0
0.8
1.5
2.1
1.4
2.7
1.5
1.1
1.2
1.7
2.2
2.1
1.6
1.2
1.6
1.8
1.1
2.3
1.2
1.6
1.4
1.7
1.6
1.0
1.5
1.8
2.0
1.5
3.0
1.7
2.0
1.7
2.4
2.7
1.5
1.3
2.2
3.3
1.4
2.1
2.0
1.7
1.4
2.1
1.7
2.3
1.2
1.7
1.5
1.7
1.7
2.3
1.8
1.5
2.7
2.3
3.0
2.4
2.4
2.9
1.7
1.5
1.0
2.5
1.7
1.9
2.0
1.8
1.4
2.1
1.6
2.3
1.7
2.3
2.7
1.8
1.4
1.6
1.7
1.2
2.5
1.5
1.9
1.4
1.9
1.5
1.9
1.4
1.7
1.5
1.5
1.6
2.1
1.7
2.3
1.0
1.5
1.5
1.4
0.9
2.8
1.6
2.1
1.8
1.7
2.3
1.8
2.0
1.3
2.1
2.0
0.5
1.2
1.2
2.2
2.2
0.8
1.2
1.8
1.0
1.9
2.0
1.7
1.9
2.5
1.1
2.2
1.1
1.4
1.4
1.7
2.1
1.4
2.0
1.9
1.7
2.5
1.2
0.9
1.2
2.2
2.9
2.5
2.0
2.1
2.0
1.8
2.0
2.1
2.0
1.5
1.5
2.7
1.8
2.6
1.4
1.9
2.6
1.5
2.1
1.6
2.2
2.0
1.5
2.1
1.8
1.9
2.0
1.8
0.9
2.0
4.6
3.6
1.6
1.4
1.3
2.0
2.9
1.3
2.3
1.7
1.5
3.1
1.8
1.4
1.7
2.9
1.9
1.2
3.0
1.7
2.5
1.3
4.2
1.4
1.6
2.2
2.2
1.5
1.6
1.7
1.7
1.3
2.1
3.1
2.6
1.6
1.5
1.7
1.0
1.0
1.9
1.4
1.2
0.9
2.2
1.6
1.4
2.2
2.2
1.2
1.5
1.2
1.9
1.3
2.1
0.9
1.2
1.3
1.5
1.5
2.4
2.4
2.3
1.9
2.0
2.2
0.8
1.8
1.9
1.1
1.2
1.5
3.4
1.6
1.5
1.4
0.9
1.7
1.5