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

SHAPE version 2b, with bootstrap Confidence Interval of Hazard Parameters

Browse Source
This commit is contained in:
Konstantinos Leptokaropoulos
2020-06-08 13:02:54 +02:00
parent 1177e49e6d
commit eab8aa62ef
41 changed files with 13059 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
SHAPE_Package/!Read_me.txt
View File

@@ -12,11 +12,15 @@ SHAPE_ver2 - is a Wrapper version (inputs are set in the wrapper script and onc
REQUIREMENTS: Matlab R2017b or later, with installed Toolboxes:
- 'Statistics and Machine Learning'
SHAPE_ver2b- is identical to SHAPE_ver2, but also estimates 95% bootstrap confidence interval
of Hazard Parameters and b-value
>> Please refer to the "READ_ME_SHAPE_ver*.pdf" file in each directory for further instructions on
step-by-step implementation of SHAPE as well as for data format specifications and general tips.
>> Your feedback is welcome! contact Dr K. Leptokaropoulos for questions, suggestions or comments:
- kleptoka@igf.edu.pl
>> Your feedback is welcome! contact Dr Kostas Leptokaropoulos for questions, suggestions or comments:
- yebeneka@gmail.com, kleptoka@igf.edu.pl
>> Please acknowledge any use of SHAPE in your work, by citing:
- Leptokaropoulos, K. and S. Lasocki (2020), SHAPE: A MATLAB Software Package for Time-dependent

8
SHAPE_Package/!What's_New!.txt Executable file
View File

@@ -0,0 +1,8 @@
>>>> 08 JUNE 2020
Please find the new version SHAPE_ver2b, which in addition to all functionalities included in SHAPE_ver2
also estimates the 95% bootstrap confidence interval of Mean Return Period, Exceedance Probability, and
(for parametric models) b-value.
The package is available in this repository together with sample data and complete user guide.

6947
SHAPE_Package/SHAPE_ver2b.0/CATALOGS/ST2_SEIS_Data.txt Executable file
View File

File diff suppressed because it is too large Load Diff

1
SHAPE_Package/SHAPE_ver2b.0/CATALOGS/ST2_SEIS_Fields.txt Executable file
View File

@@ -0,0 +1 @@
Time ML Lat Long Depth

717
SHAPE_Package/SHAPE_ver2b.0/PRODUCTION_DATA/ST2_PROD_Data.txt Executable file
View File

@@ -0,0 +1,717 @@
2013 08 23 00 00 00 139.89
2013 08 24 00 00 00 140.04
2013 08 25 00 00 00 140.32
2013 08 26 00 00 00 140.44
2013 08 27 00 00 00 140.34
2013 08 28 00 00 00 140.29
2013 08 29 00 00 00 140.42
2013 08 30 00 00 00 140.34
2013 08 31 00 00 00 140.21
2013 09 01 00 00 00 140.2
2013 09 02 00 00 00 140.3
2013 09 03 00 00 00 140.22
2013 09 04 00 00 00 140.16
2013 09 05 00 00 00 140.18
2013 09 06 00 00 00 140.19
2013 09 07 00 00 00 140.33
2013 09 08 00 00 00 140.32
2013 09 09 00 00 00 140.33
2013 09 10 00 00 00 140.38
2013 09 11 00 00 00 140.38
2013 09 12 00 00 00 140.92
2013 09 13 00 00 00 141.03
2013 09 14 00 00 00 140.95
2013 09 15 00 00 00 141.65
2013 09 16 00 00 00 141.27
2013 09 17 00 00 00 141.01
2013 09 18 00 00 00 142.7
2013 09 19 00 00 00 144.12
2013 09 20 00 00 00 144.04
2013 09 21 00 00 00 143.65
2013 09 22 00 00 00 143.16
2013 09 23 00 00 00 142.71
2013 09 24 00 00 00 142.21
2013 09 25 00 00 00 141.48
2013 09 26 00 00 00 141.46
2013 09 27 00 00 00 141.5
2013 09 28 00 00 00 141.29
2013 09 29 00 00 00 141.03
2013 09 30 00 00 00 140.56
2013 10 01 00 00 00 140.2
2013 10 02 00 00 00 145.84
2013 10 03 00 00 00 151.98
2013 10 04 00 00 00 155.71
2013 10 05 00 00 00 157.51
2013 10 06 00 00 00 158.13
2013 10 07 00 00 00 158.26
2013 10 08 00 00 00 158.1
2013 10 09 00 00 00 157.78
2013 10 10 00 00 00 157.36
2013 10 11 00 00 00 156.84
2013 10 12 00 00 00 156.28
2013 10 13 00 00 00 155.69
2013 10 14 00 00 00 154.99
2013 10 15 00 00 00 157.77
2013 10 16 00 00 00 158.99
2013 10 17 00 00 00 158.86
2013 10 18 00 00 00 158.67
2013 10 19 00 00 00 158.33
2013 10 20 00 00 00 158.13
2013 10 21 00 00 00 158.12
2013 10 22 00 00 00 157.84
2013 10 23 00 00 00 157.43
2013 10 24 00 00 00 156.91
2013 10 25 00 00 00 156.32
2013 10 26 00 00 00 155.68
2013 10 27 00 00 00 155.03
2013 10 28 00 00 00 154.33
2013 10 29 00 00 00 153.55
2013 10 30 00 00 00 152.81
2013 10 31 00 00 00 152.03
2013 11 01 00 00 00 151.16
2013 11 02 00 00 00 151.23
2013 11 03 00 00 00 149.37
2013 11 04 00 00 00 148.41
2013 11 05 00 00 00 147.61
2013 11 06 00 00 00 149.06
2013 11 07 00 00 00 154.63
2013 11 08 00 00 00 157.7
2013 11 09 00 00 00 158.46
2013 11 10 00 00 00 159.1
2013 11 11 00 00 00 159.46
2013 11 12 00 00 00 159.39
2013 11 13 00 00 00 159.14
2013 11 14 00 00 00 158.91
2013 11 15 00 00 00 163.76
2013 11 16 00 00 00 164.31
2013 11 17 00 00 00 164.34
2013 11 18 00 00 00 163.4
2013 11 19 00 00 00 162.6
2013 11 20 00 00 00 162.52
2013 11 21 00 00 00 162.12
2013 11 22 00 00 00 161.8
2013 11 23 00 00 00 161.6
2013 11 24 00 00 00 161.41
2013 11 25 00 00 00 161.21
2013 11 26 00 00 00 161.05
2013 11 27 00 00 00 160.83
2013 11 28 00 00 00 160.6
2013 11 29 00 00 00 160.97
2013 11 30 00 00 00 161.24
2013 12 01 00 00 00 161.25
2013 12 02 00 00 00 161.13
2013 12 03 00 00 00 160.96
2013 12 04 00 00 00 160.69
2013 12 05 00 00 00 160.33
2013 12 06 00 00 00 159.89
2013 12 07 00 00 00 159.37
2013 12 08 00 00 00 158.79
2013 12 09 00 00 00 158.28
2013 12 10 00 00 00 157.97
2013 12 11 00 00 00 157.95
2013 12 12 00 00 00 158.13
2013 12 13 00 00 00 158.58
2013 12 14 00 00 00 158.71
2013 12 15 00 00 00 158.99
2013 12 16 00 00 00 158.37
2013 12 17 00 00 00 159.5
2013 12 18 00 00 00 159.72
2013 12 19 00 00 00 159.95
2013 12 20 00 00 00 160.39
2013 12 21 00 00 00 160.82
2013 12 22 00 00 00 161.25
2013 12 23 00 00 00 161.54
2013 12 24 00 00 00 161.99
2013 12 25 00 00 00 162.37
2013 12 27 00 00 00 163.07
2013 12 28 00 00 00 163.41
2013 12 29 00 00 00 163.73
2013 12 30 00 00 00 164.04
2014 01 01 00 00 00 164.42
2014 01 02 00 00 00 163.74
2014 01 03 00 00 00 164.7
2014 01 04 00 00 00 164.74
2014 01 05 00 00 00 164.78
2014 01 06 00 00 00 164.9
2014 01 07 00 00 00 164.98
2014 01 08 00 00 00 164.96
2014 01 09 00 00 00 164.89
2014 01 10 00 00 00 164.88
2014 01 11 00 00 00 164.94
2014 01 13 00 00 00 165.34
2014 01 14 00 00 00 165.69
2014 01 15 00 00 00 165.82
2014 01 16 00 00 00 165.85
2014 01 17 00 00 00 165.74
2014 01 18 00 00 00 165.66
2014 01 19 00 00 00 165.58
2014 01 20 00 00 00 165.59
2014 01 21 00 00 00 165.62
2014 01 22 00 00 00 165.66
2014 01 23 00 00 00 165.73
2014 01 24 00 00 00 165.63
2014 01 25 00 00 00 165.53
2014 01 26 00 00 00 165.45
2014 01 27 00 00 00 165.38
2014 01 28 00 00 00 165.4
2014 01 29 00 00 00 165.55
2014 01 30 00 00 00 165.62
2014 01 31 00 00 00 165.58
2014 02 01 00 00 00 165.54
2014 02 02 00 00 00 165.49
2014 02 03 00 00 00 165.51
2014 02 04 00 00 00 165.43
2014 02 05 00 00 00 165.35
2014 02 06 00 00 00 165.3
2014 02 07 00 00 00 165.23
2014 02 08 00 00 00 165.07
2014 02 09 00 00 00 164.99
2014 02 10 00 00 00 164.97
2014 02 11 00 00 00 164.94
2014 02 12 00 00 00 164.95
2014 02 13 00 00 00 164.91
2014 02 14 00 00 00 164.93
2014 02 15 00 00 00 164.91
2014 02 16 00 00 00 164.7
2014 02 17 00 00 00 164.62
2014 02 18 00 00 00 164.49
2014 02 19 00 00 00 164.39
2014 02 20 00 00 00 164.38
2014 02 21 00 00 00 164.38
2014 02 22 00 00 00 164.3
2014 02 23 00 00 00 164.25
2014 02 24 00 00 00 164.24
2014 02 25 00 00 00 164.19
2014 02 26 00 00 00 164.16
2014 02 27 00 00 00 164.14
2014 02 28 00 00 00 164.05
2014 03 01 00 00 00 163.96
2014 03 02 00 00 00 163.99
2014 03 03 00 00 00 163.96
2014 03 04 00 00 00 163.93
2014 03 05 00 00 00 163.99
2014 03 06 00 00 00 164.04
2014 03 07 00 00 00 164
2014 03 08 00 00 00 163.85
2014 03 09 00 00 00 163.75
2014 03 10 00 00 00 163.66
2014 03 11 00 00 00 163.54
2014 03 12 00 00 00 163.37
2014 03 13 00 00 00 163.16
2014 03 14 00 00 00 162.94
2014 03 15 00 00 00 162.73
2014 03 16 00 00 00 162.75
2014 03 17 00 00 00 162.72
2014 03 18 00 00 00 162.26
2014 03 19 00 00 00 161.68
2014 03 20 00 00 00 161.27
2014 03 21 00 00 00 160.99
2014 03 22 00 00 00 160.67
2014 03 23 00 00 00 160.59
2014 03 24 00 00 00 160.76
2014 03 25 00 00 00 160.9
2014 03 26 00 00 00 161.01
2014 03 27 00 00 00 161.13
2014 03 28 00 00 00 161.11
2014 03 29 00 00 00 160.87
2014 03 30 00 00 00 160.62
2014 03 31 00 00 00 160.37
2014 04 01 00 00 00 160.12
2014 04 02 00 00 00 160.03
2014 04 03 00 00 00 160.15
2014 04 04 00 00 00 160.29
2014 04 05 00 00 00 160.5
2014 04 06 00 00 00 160.78
2014 04 07 00 00 00 160.86
2014 04 08 00 00 00 160.78
2014 04 09 00 00 00 160.69
2014 04 10 00 00 00 160.62
2014 04 11 00 00 00 160.36
2014 04 12 00 00 00 160.38
2014 04 13 00 00 00 160.51
2014 04 14 00 00 00 160.63
2014 04 15 00 00 00 160.73
2014 04 16 00 00 00 160.83
2014 04 17 00 00 00 160.8
2014 04 18 00 00 00 160.64
2014 04 19 00 00 00 160.44
2014 04 20 00 00 00 160.34
2014 04 21 00 00 00 160.22
2014 04 22 00 00 00 159.49
2014 04 23 00 00 00 158.81
2014 04 24 00 00 00 158.2
2014 04 25 00 00 00 158.25
2014 04 26 00 00 00 158.34
2014 04 27 00 00 00 158.45
2014 04 28 00 00 00 158.45
2014 04 29 00 00 00 158.59
2014 04 30 00 00 00 158.64
2014 05 01 00 00 00 158.74
2014 05 02 00 00 00 159.12
2014 05 03 00 00 00 159.52
2014 05 04 00 00 00 159.91
2014 05 05 00 00 00 160.2
2014 05 06 00 00 00 160.5
2014 05 07 00 00 00 160.82
2014 05 08 00 00 00 161.04
2014 05 09 00 00 00 161.27
2014 05 10 00 00 00 161.36
2014 05 11 00 00 00 161.5
2014 05 12 00 00 00 161.63
2014 05 13 00 00 00 161.83
2014 05 14 00 00 00 161.96
2014 05 15 00 00 00 161.89
2014 05 16 00 00 00 161.62
2014 05 17 00 00 00 161.59
2014 05 18 00 00 00 161.41
2014 05 19 00 00 00 161.2
2014 05 20 00 00 00 161.02
2014 05 21 00 00 00 160.81
2014 05 22 00 00 00 160.4
2014 05 23 00 00 00 159.66
2014 05 24 00 00 00 159.09
2014 05 25 00 00 00 158.79
2014 05 26 00 00 00 158.38
2014 05 27 00 00 00 157.86
2014 05 28 00 00 00 157.05
2014 05 29 00 00 00 156.29
2014 05 30 00 00 00 155.67
2014 05 31 00 00 00 155.03
2014 06 01 00 00 00 154.47
2014 06 02 00 00 00 154.21
2014 06 03 00 00 00 153.67
2014 06 04 00 00 00 153.02
2014 06 05 00 00 00 152.59
2014 06 06 00 00 00 152.23
2014 06 07 00 00 00 151.8
2014 06 08 00 00 00 151.61
2014 06 09 00 00 00 151.67
2014 06 10 00 00 00 151.38
2014 06 11 00 00 00 150.87
2014 06 12 00 00 00 150.24
2014 06 13 00 00 00 149.65
2014 06 14 00 00 00 149.27
2014 06 15 00 00 00 148.86
2014 06 16 00 00 00 148.54
2014 06 17 00 00 00 148.33
2014 06 18 00 00 00 148.48
2014 06 19 00 00 00 148.58
2014 06 20 00 00 00 148.56
2014 06 21 00 00 00 148.34
2014 06 22 00 00 00 147.96
2014 06 23 00 00 00 147.51
2014 06 24 00 00 00 147.04
2014 06 25 00 00 00 146.6
2014 06 26 00 00 00 146.17
2014 06 27 00 00 00 145.89
2014 06 28 00 00 00 145.96
2014 06 29 00 00 00 146.13
2014 06 30 00 00 00 146.39
2014 07 01 00 00 00 146.39
2014 07 02 00 00 00 146.4
2014 07 03 00 00 00 145.64
2014 07 04 00 00 00 145.25
2014 07 05 00 00 00 145.09
2014 07 06 00 00 00 145.28
2014 07 07 00 00 00 144.68
2014 07 08 00 00 00 144.73
2014 07 09 00 00 00 144.84
2014 07 10 00 00 00 145
2014 07 11 00 00 00 144.97
2014 07 12 00 00 00 144.6
2014 07 13 00 00 00 144.3
2014 07 14 00 00 00 144
2014 07 15 00 00 00 143.66
2014 07 16 00 00 00 143.24
2014 07 17 00 00 00 143.02
2014 07 18 00 00 00 143.22
2014 07 19 00 00 00 143.38
2014 07 20 00 00 00 143.47
2014 07 21 00 00 00 143.43
2014 07 22 00 00 00 143.01
2014 07 23 00 00 00 142.84
2014 07 24 00 00 00 143.36
2014 07 25 00 00 00 142.12
2014 07 26 00 00 00 141.72
2014 07 27 00 00 00 141.6
2014 07 28 00 00 00 141.95
2014 07 29 00 00 00 142.23
2014 07 30 00 00 00 142.76
2014 07 31 00 00 00 142.96
2014 08 01 00 00 00 143.06
2014 08 02 00 00 00 142.86
2014 08 03 00 00 00 142.63
2014 08 04 00 00 00 142.37
2014 08 05 00 00 00 142.08
2014 08 06 00 00 00 141.96
2014 08 07 00 00 00 142.19
2014 08 08 00 00 00 142.4
2014 08 09 00 00 00 142.58
2014 08 10 00 00 00 142.57
2014 08 11 00 00 00 142.17
2014 08 12 00 00 00 141.79
2014 08 13 00 00 00 141.43
2014 08 14 00 00 00 141.05
2014 08 15 00 00 00 140.7
2014 08 16 00 00 00 140.68
2014 08 17 00 00 00 140.93
2014 08 18 00 00 00 141.13
2014 08 19 00 00 00 141.1
2014 08 20 00 00 00 140.75
2014 08 21 00 00 00 140.69
2014 08 22 00 00 00 141.01
2014 08 23 00 00 00 141.33
2014 08 24 00 00 00 141.6
2014 08 25 00 00 00 141.8
2014 08 26 00 00 00 142.06
2014 08 27 00 00 00 142.39
2014 08 28 00 00 00 142.75
2014 08 29 00 00 00 143.02
2014 08 30 00 00 00 143.44
2014 08 31 00 00 00 143.87
2014 09 01 00 00 00 145.92
2014 09 02 00 00 00 146.24
2014 09 03 00 00 00 146.81
2014 09 04 00 00 00 147.26
2014 09 05 00 00 00 147.58
2014 09 06 00 00 00 147.71
2014 09 07 00 00 00 147.87
2014 09 08 00 00 00 148.08
2014 09 09 00 00 00 148.46
2014 09 10 00 00 00 148.58
2014 09 11 00 00 00 148.15
2014 09 12 00 00 00 147.39
2014 09 13 00 00 00 146.9
2014 09 14 00 00 00 146.85
2014 09 15 00 00 00 146.66
2014 09 16 00 00 00 145.9
2014 09 17 00 00 00 145.4
2014 09 18 00 00 00 145.19
2014 09 19 00 00 00 145.01
2014 09 20 00 00 00 145.11
2014 09 21 00 00 00 145.44
2014 09 22 00 00 00 145.72
2014 09 23 00 00 00 145.87
2014 09 24 00 00 00 146.57
2014 09 25 00 00 00 146.94
2014 09 26 00 00 00 147.25
2014 09 27 00 00 00 147.61
2014 09 28 00 00 00 147.92
2014 09 29 00 00 00 148.23
2014 09 30 00 00 00 148.53
2014 10 01 00 00 00 148.68
2014 10 02 00 00 00 148.67
2014 10 03 00 00 00 148.35
2014 10 04 00 00 00 147.9
2014 10 05 00 00 00 147.96
2014 10 06 00 00 00 148.34
2014 10 07 00 00 00 148.53
2014 10 08 00 00 00 149.56
2014 10 09 00 00 00 149.92
2014 10 10 00 00 00 149.94
2014 10 11 00 00 00 149.82
2014 10 12 00 00 00 149.55
2014 10 13 00 00 00 149.43
2014 10 14 00 00 00 148.95
2014 10 15 00 00 00 148.6
2014 10 16 00 00 00 148.33
2014 10 17 00 00 00 148.31
2014 10 18 00 00 00 148.4
2014 10 19 00 00 00 149.27
2014 10 20 00 00 00 149.57
2014 10 21 00 00 00 149.02
2014 10 22 00 00 00 148.84
2014 10 23 00 00 00 149.12
2014 10 24 00 00 00 149.12
2014 10 25 00 00 00 149.37
2014 10 26 00 00 00 150.5
2014 10 27 00 00 00 151.23
2014 10 28 00 00 00 151.17
2014 10 29 00 00 00 150.73
2014 10 30 00 00 00 150.42
2014 10 31 00 00 00 150.19
2014 11 01 00 00 00 150.07
2014 11 02 00 00 00 149.75
2014 11 03 00 00 00 149.66
2014 11 04 00 00 00 149.23
2014 11 05 00 00 00 148.96
2014 11 06 00 00 00 149.09
2014 11 07 00 00 00 149.34
2014 11 08 00 00 00 149.84
2014 11 09 00 00 00 150.31
2014 11 10 00 00 00 150.63
2014 11 11 00 00 00 150.76
2014 11 12 00 00 00 150.89
2014 11 13 00 00 00 151.55
2014 11 14 00 00 00 153.94
2014 11 15 00 00 00 156.7
2014 11 16 00 00 00 158.48
2014 11 17 00 00 00 159.22
2014 11 18 00 00 00 159.51
2014 11 19 00 00 00 159.47
2014 11 20 00 00 00 159.62
2014 11 21 00 00 00 159.73
2014 11 22 00 00 00 159.78
2014 11 23 00 00 00 160.01
2014 11 24 00 00 00 160.38
2014 11 25 00 00 00 160.44
2014 11 26 00 00 00 160.52
2014 11 27 00 00 00 160.65
2014 11 28 00 00 00 160.78
2014 11 29 00 00 00 160.98
2014 11 30 00 00 00 161.45
2014 12 01 00 00 00 162.55
2014 12 02 00 00 00 162.91
2014 12 03 00 00 00 163.15
2014 12 04 00 00 00 163.01
2014 12 05 00 00 00 163.36
2014 12 06 00 00 00 163.98
2014 12 07 00 00 00 164.71
2014 12 08 00 00 00 164.92
2014 12 09 00 00 00 165.15
2014 12 10 00 00 00 165.02
2014 12 11 00 00 00 164.78
2014 12 12 00 00 00 165.08
2014 12 13 00 00 00 165.33
2014 12 14 00 00 00 165.39
2014 12 15 00 00 00 165.25
2014 12 16 00 00 00 165.36
2014 12 17 00 00 00 165.8
2014 12 18 00 00 00 165.82
2014 12 19 00 00 00 165.67
2014 12 20 00 00 00 165.41
2014 12 21 00 00 00 165.2
2014 12 22 00 00 00 165.32
2014 12 23 00 00 00 165.57
2014 12 24 00 00 00 165.36
2014 12 25 00 00 00 165.31
2014 12 26 00 00 00 165.42
2014 12 27 00 00 00 165.75
2014 12 28 00 00 00 165.99
2014 12 29 00 00 00 166
2014 12 30 00 00 00 166.06
2014 12 31 00 00 00 166.79
2015 01 01 00 00 00 167.54
2015 01 02 00 00 00 168.16
2015 01 03 00 00 00 168.72
2015 01 04 00 00 00 169.3
2015 01 05 00 00 00 169.77
2015 01 06 00 00 00 170.17
2015 01 07 00 00 00 170.74
2015 01 08 00 00 00 171.51
2015 01 09 00 00 00 171.91
2015 01 10 00 00 00 172
2015 01 11 00 00 00 172
2015 01 12 00 00 00 171.97
2015 01 13 00 00 00 171.77
2015 01 14 00 00 00 171.53
2015 01 15 00 00 00 171.22
2015 01 16 00 00 00 170.97
2015 01 17 00 00 00 170.68
2015 01 18 00 00 00 170.65
2015 01 19 00 00 00 170.78
2015 01 20 00 00 00 170.93
2015 01 21 00 00 00 171.05
2015 01 22 00 00 00 171.07
2015 01 23 00 00 00 171.04
2015 01 24 00 00 00 170.99
2015 01 25 00 00 00 170.94
2015 01 26 00 00 00 170.91
2015 01 27 00 00 00 170.95
2015 01 28 00 00 00 170.93
2015 01 29 00 00 00 171.02
2015 01 30 00 00 00 171.05
2015 01 31 00 00 00 171.1
2015 02 01 00 00 00 171.28
2015 02 02 00 00 00 171.57
2015 02 03 00 00 00 171.29
2015 02 04 00 00 00 170.9
2015 02 05 00 00 00 170.82
2015 02 06 00 00 00 170.77
2015 02 07 00 00 00 170.73
2015 02 08 00 00 00 170.66
2015 02 09 00 00 00 170.6
2015 02 10 00 00 00 170.53
2015 02 11 00 00 00 170.4
2015 02 12 00 00 00 170.37
2015 02 13 00 00 00 170.34
2015 02 14 00 00 00 170.4
2015 02 15 00 00 00 170.58
2015 02 16 00 00 00 170.63
2015 02 17 00 00 00 170.64
2015 02 18 00 00 00 170.8
2015 02 19 00 00 00 170.88
2015 02 20 00 00 00 171.11
2015 02 21 00 00 00 171.33
2015 02 22 00 00 00 171.59
2015 02 23 00 00 00 171.72
2015 02 24 00 00 00 171.69
2015 02 25 00 00 00 171.53
2015 02 26 00 00 00 171.47
2015 02 27 00 00 00 171.46
2015 02 28 00 00 00 171.38
2015 03 01 00 00 00 171.37
2015 03 02 00 00 00 171.35
2015 03 03 00 00 00 171.27
2015 03 04 00 00 00 171.15
2015 03 05 00 00 00 171.04
2015 03 06 00 00 00 170.9
2015 03 07 00 00 00 170.82
2015 03 08 00 00 00 170.7
2015 03 09 00 00 00 170.57
2015 03 10 00 00 00 170.55
2015 03 11 00 00 00 170.69
2015 03 12 00 00 00 170.86
2015 03 13 00 00 00 170.95
2015 03 14 00 00 00 170.88
2015 03 15 00 00 00 170.76
2015 03 16 00 00 00 170.5
2015 03 17 00 00 00 170.29
2015 03 18 00 00 00 170.17
2015 03 19 00 00 00 169.91
2015 03 20 00 00 00 169.5
2015 03 21 00 00 00 169.28
2015 04 01 00 00 00 170.68
2015 04 02 00 00 00 170.36
2015 04 03 00 00 00 170.02
2015 04 04 00 00 00 169.7
2015 04 05 00 00 00 169.55
2015 04 06 00 00 00 169.41
2015 04 07 00 00 00 169.04
2015 04 08 00 00 00 168.37
2015 04 09 00 00 00 167.64
2015 04 10 00 00 00 167.32
2015 04 11 00 00 00 167.49
2015 04 12 00 00 00 167.63
2015 04 13 00 00 00 167.87
2015 04 14 00 00 00 168.05
2015 04 15 00 00 00 168.21
2015 04 16 00 00 00 168.34
2015 04 17 00 00 00 168.48
2015 04 18 00 00 00 168.61
2015 04 19 00 00 00 168.75
2015 04 20 00 00 00 168.89
2015 04 21 00 00 00 169.01
2015 04 22 00 00 00 169.18
2015 04 23 00 00 00 169.09
2015 04 24 00 00 00 168.78
2015 04 25 00 00 00 168.29
2015 04 26 00 00 00 167.85
2015 04 27 00 00 00 167.34
2015 04 28 00 00 00 167.01
2015 04 29 00 00 00 166.79
2015 04 30 00 00 00 166.7
2015 05 01 00 00 00 166.83
2015 05 02 00 00 00 166.65
2015 05 03 00 00 00 166.01
2015 05 04 00 00 00 165.45
2015 05 05 00 00 00 164.61
2015 05 06 00 00 00 163.86
2015 05 07 00 00 00 163.36
2015 05 08 00 00 00 163.02
2015 05 09 00 00 00 162.73
2015 05 10 00 00 00 162.81
2015 05 11 00 00 00 162.93
2015 05 12 00 00 00 162.76
2015 05 13 00 00 00 162.54
2015 05 14 00 00 00 162.36
2015 05 15 00 00 00 162.19
2015 05 16 00 00 00 162.03
2015 05 17 00 00 00 161.93
2015 05 18 00 00 00 161.79
2015 05 19 00 00 00 161.38
2015 05 20 00 00 00 160.79
2015 05 21 00 00 00 160.15
2015 05 22 00 00 00 159.43
2015 05 23 00 00 00 158.84
2015 05 24 00 00 00 158.38
2015 05 25 00 00 00 158
2015 05 26 00 00 00 157.48
2015 05 27 00 00 00 157.09
2015 05 28 00 00 00 156.77
2015 05 29 00 00 00 156.41
2015 05 30 00 00 00 156.06
2015 05 31 00 00 00 155.75
2015 06 01 00 00 00 155.45
2015 06 02 00 00 00 155.02
2015 06 03 00 00 00 154.79
2015 06 04 00 00 00 154.45
2015 06 05 00 00 00 154.18
2015 06 06 00 00 00 153.97
2015 06 07 00 00 00 153.78
2015 06 08 00 00 00 153.53
2015 06 09 00 00 00 153.28
2015 06 10 00 00 00 153.1
2015 06 11 00 00 00 152.48
2015 06 12 00 00 00 152.02
2015 06 13 00 00 00 152.06
2015 06 14 00 00 00 151.98
2015 06 15 00 00 00 152.06
2015 06 16 00 00 00 151.69
2015 06 17 00 00 00 151.11
2015 06 18 00 00 00 150.89
2015 06 19 00 00 00 150.73
2015 06 20 00 00 00 150.5
2015 06 21 00 00 00 150.41
2015 06 22 00 00 00 150.56
2015 06 23 00 00 00 150.56
2015 06 24 00 00 00 150.44
2015 06 25 00 00 00 150.18
2015 06 26 00 00 00 149.74
2015 06 28 00 00 00 149.24
2015 06 29 00 00 00 148.92
2015 06 30 00 00 00 148.5
2015 07 01 00 00 00 148.19
2015 07 02 00 00 00 148.16
2015 07 03 00 00 00 147.92
2015 07 04 00 00 00 147.4
2015 07 05 00 00 00 147.07
2015 07 06 00 00 00 146.79
2015 07 07 00 00 00 146.46
2015 07 08 00 00 00 146.11
2015 07 09 00 00 00 145.9
2015 07 10 00 00 00 146.13
2015 07 11 00 00 00 146.38
2015 07 12 00 00 00 146.58
2015 07 13 00 00 00 146.57
2015 07 14 00 00 00 146.36
2015 07 15 00 00 00 146.08
2015 07 16 00 00 00 145.89
2015 07 17 00 00 00 146.06
2015 07 18 00 00 00 146.2
2015 07 19 00 00 00 146.45
2015 07 21 00 00 00 146.53
2015 07 22 00 00 00 146.16
2015 07 23 00 00 00 145.81
2015 07 24 00 00 00 145.4
2015 07 25 00 00 00 145.02
2015 07 26 00 00 00 144.61
2015 07 27 00 00 00 144.41
2015 07 28 00 00 00 144.51
2015 07 29 00 00 00 144.65
2015 07 30 00 00 00 144.82
2015 07 31 00 00 00 144.75
2015 08 01 00 00 00 144.38
2015 08 02 00 00 00 144.01
2015 08 03 00 00 00 143.67
2015 08 04 00 00 00 143.49
2015 08 05 00 00 00 143.19
2015 08 06 00 00 00 143.28
2015 08 07 00 00 00 144.03
2015 08 08 00 00 00 144.44
2015 08 09 00 00 00 144.81
2015 08 10 00 00 00 145.19
2015 08 11 00 00 00 145.23
2015 08 12 00 00 00 145.14
2015 08 13 00 00 00 144.96
2015 08 14 00 00 00 144.69
2015 08 15 00 00 00 144.43
2015 08 16 00 00 00 144.29
2015 08 17 00 00 00 144.5
2015 08 18 00 00 00 144.68
2015 08 19 00 00 00 144.82
2015 08 20 00 00 00 144.8
2015 08 21 00 00 00 144.59
2015 08 22 00 00 00 144.34
2015 08 23 00 00 00 144.07
2015 08 24 00 00 00 143.76

1
SHAPE_Package/SHAPE_ver2b.0/PRODUCTION_DATA/ST2_PROD_Fields.txt Executable file
View File

@@ -0,0 +1 @@
Date Water_Level

BIN
SHAPE_Package/SHAPE_ver2b.0/READ_ME_SHAPE_ver2b.pdf Executable file
View File

Binary file not shown.

529
SHAPE_Package/SHAPE_ver2b.0/SHAPE_ver2b_1.m Executable file
View File

@@ -0,0 +1,529 @@
% PROGRAM : SHAPE [Seismic HAzard Parameters Evaluation]
% VERSION : V_2b.1 [Wrapper (fast) Standalone Version]
% LAST UPDATED: June 2020
% COMPATIBLE : Matlab version 2017b or later
% REQUIREMENTS: "Statistics and Machine Learning" Matlab Toolbox installed
% TOOLBOX : "Hazard Analysis Toolbox" developed within SERA Project
% USER GUIDE : "READ_ME_SHAPE_ver2b.pdf"
% CITE AS : Leptokaropoulos K. and S. Lasocki (2020), Seismol. Res. Lett., doi: 10.1785/0220190319
% CONTACT : kleptoka@igf.edu.pl (Dr. Kostas Leptokaropoulos)
% --------------------------------------------------------------------------------------------------------------------
% Time-and-Technology Dependent Seismic Hazard Assessment (SHA)
% --------------------------------------------------------------------------------------------------------------------
% INPUT:
% !!! ---------------------------- INPUT DATA REQUIREMENTS ----------------------------- !!!
% the program works with ASCII input data files (e.g. *.txt). The files needed are:
% > File with the parameters of seismic data [mandatory]
% > File with the parameters of production data [optional]
% > File specifying time windows for SHA analysis [optional]
% > File with the fields description of the corresponding parameters in the seismic data file
% > File with the fields description of the corresponding parameters in the production data files
% FOR DETAILS on data requirements please refer to the document:
% "READ_ME_SHAPE_ver2b.pdf"
% --------------------------------------------------------------------------------------------------------------------
% OVERVIEW: THE PROGRAM takes as input a Seismic and optionally, a Production data parameter file to provide Seismic
% Hazard Assessment for specified time-windows, following User's specifications.
% --------------------------------------------------------------------------------------------------------------------
% AUTHORS: K. Leptokaropoulos and S. Lasocki
% Based on magnitude distribution/stationary hazard parameters estimation functions originally developed by S. Lasocki
% Last Updated: 06/2020, within SERA PROJECT, EU Horizon 2020 R&I programme under grant agreement No.730900
% CURRENT VERSION: v2b.1 **** [WRAPPER (fast) STANDALONE VERSION!!]
%% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PLEASE refer to the accompanying document:
% "READ_ME_SHAPE_ver2b.pdf"
% for description of the Application and its requirements.
%% -----------------------------------------------------------------------------------------------------------------------
% DESCRIPTION: The Application performs time-dependent Seismic Hazard Analysis (SHA),
% taking into account the activity rate and the magnitude distribution of seismicity
% for selected time windows. The hazard parameters estimated are:
% 1) The Mean Return Period (MRP) of a given magnitude, M, which is defined as the
% average elapsed time between the occurrence of consecutive events of M and
% 2) The Exceedance Probability (EPR) of a given magnitude, M, within a given time
% period of length, T, which is defined as the probability of an earthquake of
% M to occur during T.
% These hazard parameters together with the 95% confidence interval are estimated for
% different time windows which are constructe upon User<EFBFBD>s particular specifications.
% 4 different magnitude distribution models can be chosen. The input files must be in
% ASCII format (e.g. *.txt). A brief description of the preparation process is given here:
% !!! THE USER SETS THE PARAMETERS IN LINES 125-143 !!! then the following steps are executed
% (please see in the script for more comments and details for each STEP)
% STEP_1. Upload Data
% STEP_2. Select Magnidue Scale
% STEP_3. Filter Data for Mc
% STEP_4. TIME WINDOWS GENERATION
% STEP_5. SHA PARAMETERS ESTIMATION
% STEP_6. Plotting (optional)
% STEP_7. save OUTPUTS
% ---------------------------------------------------------------------------------------------------------------------
%% INPUT: All input data are sufficiently explained in the script as well as while running the code
% (interaction with the user). NOTE that all input files (seismic catalog, production data,
% time windows) must be in ASCII format (i.e. *.txt). Please refer to the APPLICATION DOCUMENTATION
% for further instructions and input data requirement specifications: "READ_ME_SHAPE_ver2b.pdf"
% ----------------------------------------------------------------------------------------------------------------------
%% OUTPUT:
% <> Output Report with summary of the Results as well as data and parameters used
% <> Output Figure with the results in *.mat and *.jpeg formats (optional)
% <> Output Matlab Structure with input parameter values and output results, having
% as many cells as the number of time windows generated.
% Structure fields are:
% - Time : vector with origin times of the events included in each time window
% - M : vector with events magnitudes
% - Mmin : Completeness magnitude
% - eps : Magnitude round-off interval
% - lambd : mean activity rate
% - lambd_err : events number sufficiency (0-all parametes estiamated, 1-all parameters set as NaNs)
% - unit : Time Unit
% - method : Magnitude Distribution Model
% - b : b-value of GR law
% [applies only when "method" is set to 'GRU' or 'GRT']
% - b_0025 : these two parameters define the b-value 95% bootstrap Confidence interval
% - b_0075 : -//-
% - h : Kernel smoothing factor
% [applies only when "method" is set to 'NPU' or 'NPT']
% - xx : Background sample for kernel magnitude estimate
% [applies only when "method" is set to 'NPU' or 'NPT']
% - ambd : weigthing factors for the adaptive kernel
% [applies only when "method" is set to 'NPU' or 'NPT']
% - ierr : h convergence indicator (0-converges,1-multiple zeros, 2-no zeros)
% [applies only when "method" is set to 'NPU' or 'NPT']
% - Mmax : Upper limit of magnitude distribution (truncated)
% [applies only when "method" is set to 'GRT' or 'NPT']
% - err : Mmax convergence indicator (0-converge, 1-no converge)
% [applies only when "method" is set to 'GRT' or 'NPT']
% - PDF : 2 columns, first representing magnitudes and second the Probability Density Function of those magnitudes
% - CDF : 2 columns, first representing magnitudes and second the Cumulative Distribution Function of those magnitudes
% - MRP : Mean Return Period
% - MRP_0.025 : these two parameters define the MRP 95% bootstrap Confidence interval
% - MRP_0.975 : -//-
% - EP : Exceedance Probability
% - EP_0.025 : these two parameters define the EP 95% bootstrap Confidence interval
% - EP_0.975 : -//-
% -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
% REFERENCES:
% Kijko A, Lasocki S, Graham G (2001), Pure Appl. Geophys. 158:1655<EFBFBD>1676
% Kijko A, Sellevoll MA (1989), Bull Seismol. Soc. Am. 79:645<EFBFBD>654
% Lasocki S (2017), Chapter 11.3 in Rockburst Mechanisms, Oxford, pp 366-380
% Lasocki S, Urban P (2011), Acta Geophys. 59:659<EFBFBD>673
% Lasocki S, Orlecka-Sikora B (2008), Tectonophysics 456:28<EFBFBD>37
% Leptokaropoulos K, Staszek M, Cielesta S, Urban P, Olszewska D, Lizurek G (2017), Acta Geophys. 65:493-505
% Leptokaropoulos K, Lasocki S (2020), Seismol. Res. Lett., doi: 10.1785/0220190319.
% ---------------------------------------------------------------------------------------------------------------------
% 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.
% -------------------------------------------------------------------------------------------------------------------
clc;clear;
close all;mkdir Outputs_SHA
% PLEASE SET INPUT ARGUMENTS [LINES 125-143]
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SEIS_DATA='ST2_SEIS_Data.txt'; % Seismic Data File - NOTE: SEIS_DATA=[] is valid as well
SEIS_FIELDS='ST2_SEIS_Fields.txt'; % Seismic Data Fields File
PROD_DATA='ST2_PROD_Data.txt'; % Production Data (non-seismic) File - NOTE: PROD_DATA=[] is valid as well
PROD_FIELDS='ST2_PROD_Fields.txt'; % Production Data (non-seismic) Fields File
PROD_FIELD=2; % field (column) corresponding to a selected Production Parameter
MScale='ML'; % Magnitude Scale (e.g. 'ML', 'Mw' etc). NOTE: MScale=[] is also valid
Mc=1.0; % Select Mc
Mmax=[]; % Valid for GRT and NPT/ if Mmax=[], it is calculated internally
Nsynth=1000; % number of trials to estimated Mmax bias/ if Nsynth=[], no bias is estimated
winmode='Time'; % Select MODE for windows creation: 'Time' or 'Events' or 'File'
file_n='ST2_test_timewindows.txt'; % Select file name with starting and ending time of time windows, applicable only for winmode='File'
window_size=50; % time window span (days or events, depending on "winmode")
dt=50; % time step (days)
method='GRT'; % Select M distribution model among 'GRU','GRT','NPU','NPT'
Tunit='day'; % Select time unit among 'day', 'month', 'year'
MaG=3.0; % set target Magnitude for EPR and MRP calculation
Plength=1; % set target time Period (days) for EPR calculation
Nbst=100; % number of bootstrap iterations to determine Hazard Parameters confidence interval
Plotopt='ON'; % To enable ('ON') or disable ('OFF') plotting
%% -------------------------------------------------------------------------------------------------------------------------
if method=='GRT' & isempty(Mmax)==0 & MaG>Mmax | method=='NPT' & isempty(Mmax)==0 & MaG>Mmax
error('Mmax is smaller than the target magnitude (MaG), please change input')
end
%% STEP 1: Load data
[Catalog,PROD_Data,s1]=Data_Hand_A2M_ver2...
(SEIS_DATA,SEIS_FIELDS,PROD_DATA,PROD_FIELDS,PROD_FIELD);
%% STEP 2: Select Magnitude Scale
[Ctime,Cmag]=Select_Magnitude_Scale_ver2(Catalog,MScale);
%% STEP 3: Filter data for Mc
[Ctime,Cmag,Catalog]=FiltMc_ver2(Ctime,Cmag,Catalog,s1,Mc);
%% STEP 3b: Estimate Mmax for Truncated Distributions (GRT and NPT)
cd SSH/
if strcmp(method,'GRT')==1 & isempty(Mmax)==1;
[p1,p2,p3,p4,p5,pb,Mmax,err,BIAS,SD]=TruncGR_O(Ctime,Cmag,0,Mc,Mmax,Nsynth);
Mmax=Mmax+BIAS;
elseif strcmp(method,'NPT')==1 & isempty(Mmax)==1;
[p1,p2,p3,p4,p5,p6,p7,p8,p9,Mmax,err,BIAS,SD]=Nonpar_tr_O(Ctime,Cmag,0,Mc,Mmax,Nsynth);
Mmax=Mmax+BIAS;
end
cd ../
%% STEP 4: Create Time Windows
time_windows=struct;time_windows.Time=[];
to=Ctime-Ctime(1);tmin=min(to);tmax=max(to);
switch winmode
%% TIME
case 'Time'
if window_size>tmax;n=1;warning('time window is set larger than data time span');end
n=ceil((tmax-window_size)/dt);
for i=1:n
time_windows(i).Time=Ctime(to>=(i-1)*dt & to<(i-1)*dt+window_size);
time_windows(i).M=Cmag(to>=(i-1)*dt & to<(i-1)*dt+window_size);
time_windows(i).Tstart=Ctime(1)+(i-1)*dt;
time_windows(i).Tend=Ctime(1)+(i-1)*dt+window_size;
end
%% EVENTS
case 'Events'
if window_size>numel(to);window_size=numel(to);warning('events window is set larger than given data');end
n=ceil(tmax/dt);
for i=1:n
To=find(to>=(i-1)*dt);To=To(1);
if To<=length(Ctime)-window_size+1;
time_windows(i).Time=Ctime(To:To+window_size-1);
time_windows(i).M=Cmag(To:To+window_size-1);
time_windows(i).Tstart=Ctime(1)+(i-1)*dt;
time_windows(i).Tend=max(time_windows(i).Time);
end
end
%% FILE
case 'File'
cd TIME_WINDOWS
Twindows=dlmread(file_n);
T1=Twindows(:,1);T2=Twindows(:,2);n=numel(T1);
for i=1:n
time_windows(i).Time=Ctime(Ctime>=T1(i) & Ctime<T2(i));
time_windows(i).M=Cmag(Ctime>=T1(i) & Ctime<T2(i));
time_windows(i).Tstart=T1(i);
time_windows(i).Tend=T2(i);
end
cd ../
end
%% STEP 5: ESTIMATE HAZARD PARAMETERS
if strcmp(Tunit,'day');iop=0;elseif strcmp(Tunit,'month');iop=1;elseif strcmp(Tunit,'year');iop=2;
else; error('Please set "Tunit" Parameter as "day", "month" or "year"');end
%% RUN MAGDIST
[HP] = TDHMagDistWrapper(method, time_windows, Mc, iop,Mmax,Nbst);
%% Harzard Parameters Estimate
%[MRPer,ExPr,MRPer_low,MRPer_high,ExPr_low,ExPr_high]=TDHRetPeriodExcProbWrapper(method,MaG,Plength,Mc,HP,Plength)
for i=1:length(HP)
if isnan(HP(i).CDF)
lambda(i)=HP(i).lamb;MRPer(i)=NaN;ExPr(i)=NaN;MRPer_low(i)=NaN;MRPer_high(i)=NaN;ExPr_low(i)=NaN;ExPr_high(i)=NaN;
else
[MRPer1,ExPr1,MRPer_low1,MRPer_high1,ExPr_low1,ExPr_high1]=RetPeriodExcProbWrapper(HP(i).lamb,HP(i).CDF,Plength,MaG);
lambda(i)=HP(i).lamb;
MRPer(i)=MRPer1;ExPr(i)=ExPr1;
MRPer_low(i)=MRPer_low1;MRPer_high(i)=MRPer_high1;
ExPr_low(i)=ExPr_low1;ExPr_high(i)=ExPr_high1;
end
end
%%
%% STEP 7: Ploting (optional)
Zplo_ver2b_1
%% STEP 6: Save outputs
Zsave_output_ver2b_1
%% -*-*-*-*-*-*-*-*-*-*-*-*-*- F U N C T I O N S -*-*-*-*-*-*-*-*-*-*-*-*-*-
%% ------------------------- DATA HANDLING FUNCTION --------------------------
function [Catalog,PROD_data,ss1]=Data_Hand_A2M_ver2...
(SEIS_DATA,SEIS_FIELDS,PROD_DATA,PROD_FIELDS,PROD_FIELD)
if isempty(PROD_DATA)
%% SEISMIC DATA
cd CATALOGS
SData=load(SEIS_DATA);
SFields=fileread(SEIS_FIELDS);
cd ../
Na=length(SFields);if SFields(Na)~=' ';SFields(Na+1)=' ';end
[cou,c,Datime,Catalog]=Fields_dat(SData,SFields,1);
ss1=1:length(Catalog);% ss1=SetParams(Catalog);%Catalog=Catalog(ss1);
PROD_data=[];ss=[];s2=[];ss2=[];dstr2=[];
else
%% ------------------------------------------------------------------------------------------------------------------------------------
% BOTH -SEISMIC and PRODUCTION DATA
cd CATALOGS %Seismic Data
SData=load(SEIS_DATA);
SFields=fileread(SEIS_FIELDS);
cd ../
Na=length(SFields);if SFields(Na)~=' ';SFields(Na+1)=' ';end
[cou,c,Datime,Catalog]=Fields_dat(SData,SFields,1);
ss1=1:length(Catalog);
% ss1=SetParams(Catalog);%Catalog=Catalog(ss1);
cd PRODUCTION_DATA % Production (non-Seismic) Data
OData=load(PROD_DATA);
OFields=fileread(PROD_FIELDS);
cd ../
Na1=length(OFields);if OFields(Na1)~=' ';OFields(Na1+1)=' ';end
[cou1,c1,Datime1,PROD_data]=Fields_dat(OData,OFields,2);
end
% save('Catalog_ST2_Test','Catalog')
% save('Catalog_ST2_Test','Catalog')
%% ----------------------------------------- F U N C T I O N S -----------------------------------------
function [cou,c,Datime,OUT]=Fields_dat(indata,infields,iop)
Datime=datenum(indata(:,1),indata(:,2),indata(:,3),indata(:,4),indata(:,5),indata(:,6)); % Convert time to matlab format
% Define Fields
c=1;
for i=1:length(infields)-1
if strcmp(infields(i),' ')==1;cou(i)=0;
else cou(i)=c;end
if strcmp(infields(i),' ')==0 & strcmp(infields(i+1),' ')==1;c=c+1;end
end
if iop==1
OUT(1).field='Occurrence_Time';OUT(1).val=Datime;
elseif iop==2
OUT(1).field='Production_Time';OUT(1).val=Datime;
end
for i=2:c-1
OUT(i).field=infields(cou==i);
OUT(i).val=indata(:,i+5);
end
%Set Field Type for Magnitude Recognition
for i=1:size(OUT,2)
if strcmp(OUT(i).field,'ML') || strcmp(OUT(i).field,'Mw') || strcmp(OUT(i).field,'M') ...
|| strcmp(OUT(i).field,'Ms') || strcmp(OUT(i).field,'mb') || strcmp(OUT(i).field,'Md')
OUT(i).fieldType='Magnitude';
else
OUT(i).fieldType=[];
end
end
end
%% -----------------------------------------------------------------------------------------------------------
% ----------------------------------------------------------------------------------------------------------------
end
%%
function [Ctime,Cmag,Data]=FiltMc_ver2(Ctime,Cmag,Catalog,s1,Mc)
clc
id_time=findfield(Catalog,'Occurrence_Time');
%opts.Interpreter='tex';opts.Default='Yes';
%quest='Do you wish to filter Data for Mc?';
%answer=questdlg(quest,'Data Completenes','Yes','No',opts);
%if strcmp(answer,'Yes')
%% THIS HAS MOVED TO A SEPARATE FUNCTION IN THE BEGINNING OF THE APPLICATION
% cou=1;
% for i=1:length(Catalog)
% if strcmp(Catalog(i).fieldType,'Magnitude')==1
% C(cou).field=Catalog(i).field;cou=cou+1;
% end
% end
%
% % Check for no magnitude
% if cou==1;error('MyComponent:incorrectType',...
% 'No magnitude column detected!! Please check: Magnitude fields must be noted as one of the following:\nM , Mw , ML , Ms , mb , Md\n or select the entire sample for analysis');end
%
% %Select Parameters from Seismic Catalog -
% [ss1,ok]=listdlg('PromptString','Please Select M scale:',...
% 'ListString',{C.field}, 'SelectionMode','single');
%
% id=findfield(Catalog,C(ss1).field);
% Mtype=Catalog(id).field;
% Ctime=Catalog(id_time).val;
% id_M=findfield(Catalog,Mtype);
% Cmag=Catalog(id_M).val;
%%
cou=1;
for i=1:length(s1)
x=Catalog(s1(i)).val;x=x(Cmag>=Mc);
index=isnan(x)==0;
x=x(isnan(x)==0);
Data(cou).field=Catalog(s1(i)).field;
Data(cou).fieldType=Catalog(s1(i)).fieldType;
Data(cou).val=nan(size(index)); Data(cou).val(index)=x;
cou=cou+1;
end
Ctime=Ctime(Cmag>=Mc);Cmag=Cmag(Cmag>=Mc);
n=length(Ctime);
disp(['number of events: ',num2str(n)])
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 [Ctime,Cmag]=Select_Magnitude_Scale_ver2(Catalog,MType)
cou=1;
for i=1:length(Catalog)
if strcmp(Catalog(i).fieldType,'Magnitude')==1
C(cou).field=Catalog(i).field;cou=cou+1;
end
end
% Check for no magnitude
if cou==1;error('MyComponent:incorrectType',...
'No magnitude column detected!! Please check: Magnitude fields must be noted as one of the following:\nM , Mw , ML , Ms , mb , Md\n or select the entire sample for analysis');end
id=findfield(Catalog,MType);
id_time=findfield(Catalog,'Occurrence_Time');
Cmag=Catalog(id).val;Ctime=Catalog(id_time).val;
end
%% -----------!!!!!!!!!!! HAZARD PARAMETERS ESTIMATE FUNCTIONS !!!!!!!!!!!-----------
function [HP] = TDHMagDistWrapper(method, time_win_data, mmin, iop,Mmax,Nbst)
cd SSH
for j=1:size(time_win_data,2);if isempty(time_win_data(j).Time);M1(j)=NaN;else;M1(j)=max(time_win_data(j).M);end;end
Mmaxcat=max(M1);Mmax=max([Mmaxcat Mmax])
for i=1:size(time_win_data,2)
mags_vec = time_win_data(i).M;
time_vec = time_win_data(i).Time;
HP(i).mmin = mmin;
HP(i).iop = iop;
HP(i).method = method;
switch method
case 'GRU'
try
[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).b, HP(i).bCI]=UnlimitGR(time_vec, mags_vec, iop, mmin,Nbst);
HP(i).Mmax=NaN;eps=HP(i).eps;
[m, PDF_GRU, CDF_GRU]=dist_GRU_CI(mmin,Mmax+eps,eps,mmin,eps,HP(i).b,HP(i).bCI); %K29JAN2020 - magnitude PDF/CDF
HP(i).PDF=[m PDF_GRU];HP(i).CDF=[m CDF_GRU]; %K29JAN2020 - magnitude PDF/CDF
catch err
HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).b=NaN;HP(i).Mmax=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
HP(i).bCI(1)=NaN;HP(i).bCI(2)=NaN;warning('%s: %s', err.identifier, err.message);
end
case 'GRT'
try
[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).b, HP(i).Mmax, HP(i).err]=TruncGR_O(time_vec, mags_vec, iop, mmin,Mmax);
eps=HP(i).eps;[m, PDF_GRT, CDF_GRT]=dist_GRT(mmin,Mmax+eps,eps,mmin,eps,HP(i).b,Mmax); %K29JAN2020 - magnitude PDF/CDF
funct=@(mags_vec)TruncGR_O_CI(time_vec,mags_vec,iop,mmin,Mmax,[]); % K04JUN2020 GRT
HP(i).bCI=bootci(Nbst,{funct,mags_vec},'alpha',0.05); % K04JUN2020 GRT
[m, PDF_GRT, CDF_GRT]=dist_GRT_CI(mmin,Mmax+eps,eps,mmin,eps,HP(i).b,Mmax,HP(i).bCI); % K04JUN2020 GRT
HP(i).PDF=[m PDF_GRT];HP(i).CDF=[m CDF_GRT]; %K29JAN2020 - magnitude PDF/CDF
catch err
HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).b=NaN; HP(i).Mmax=NaN; HP(i).err=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
HP(i).bCI(1)=NaN;HP(i).bCI(2)=NaN;warning('%s: %s', err.identifier, err.message);
end
case 'NPU'
try
[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).ierr, HP(i).h, HP(i).xx, HP(i).ambd]=Nonpar_O(time_vec, mags_vec, iop, mmin);
HP(i).Mmax=NaN;eps=HP(i).eps;
[m, PDF_NPU, CDF_NPU]=dist_NPU(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,HP(i).xx,HP(i).ambd); %K29JAN2020 - magnitude PDF/CDF
HP(i).PDF=[m PDF_NPU];HP(i).CDF=[m CDF_NPU]; %K29JAN2020 - magnitude PDF/CDF
data=[HP(i).xx,HP(i).ambd']; % K04JUN2020 NPU
[CDF_NP]=dist_NPU_CI(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,data); % K04JUN2020 NPU
funct=@(data)dist_NPU_CI(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,data); % K04JUN2020 NPU
[CDF_CI]=bootci(Nbst,{funct,data},'alpha',0.05)'; % K04JUN2020 NPU
HP(i).CDF=[HP(i).CDF CDF_CI]; % K04JUN2020 NPU
catch err
HP(i).lamb_all=NaN; HP(i).lamb=NaN;HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).ierr=NaN; HP(i).h=NaN; HP(i).xx=[]; HP(i).ambd=[];HP(i).Mmax=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
warning('%s: %s', err.identifier, err.message);
end
case 'NPT'
try
[HP(i).lamb_all, HP(i).lamb, HP(i).lamb_err, HP(i).unit, HP(i).eps, HP(i).ierr, HP(i).h, HP(i).xx, HP(i).ambd, HP(i).Mmax, HP(i).err]=Nonpar_tr_O(time_vec, mags_vec, iop, mmin,Mmax);
eps=HP(i).eps;
[m, PDF_NPT, CDF_NPT]=dist_NPT(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,HP(i).xx,HP(i).ambd,Mmax); %K29JAN2020 - magnitude PDF/CDF
HP(i).PDF=[m PDF_NPT];HP(i).CDF=[m CDF_NPT]; %K29JAN2020 - magnitude PDF/CDF
data=[HP(i).xx,HP(i).ambd']; % K04JUN2020 NPT
[CDF_NP]=dist_NPT_CI(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,data,Mmax); % K04JUN2020 NPT
funct=@(data)dist_NPT_CI(mmin,Mmax+eps,eps,mmin,eps,HP(i).h,data,Mmax); % K04JUN2020 NPT
[CDF_CI]=bootci(Nbst,{funct,data},'alpha',0.05)'; % K04JUN2020 NPT
HP(i).CDF=[HP(i).CDF CDF_CI]; % K04JUN2020 NPU
catch err
HP(i).lamb_all=NaN; HP(i).lamb=NaN; HP(i).lamb_err=2; HP(i).unit=''; HP(i).eps=NaN; HP(i).ierr=NaN; HP(i).h=NaN; HP(i).xx=[]; HP(i).ambd=[]; HP(i).Mmax=NaN; HP(i).err=NaN;HP(i).PDF=NaN;HP(i).CDF=NaN;
warning('%s: %s', err.identifier, err.message);
end
end
% K12NOV2015
% Calculate lamb and lamb_all in case of 0, 1, or 2 events.
% It may be generalized in all cases. Rate is now calculated
% by division of the event number by the duration of the set
% - not the time difference between the first and last events.
if numel(mags_vec(mags_vec>=mmin))<3
HP(i).lamb_all=NaN; % numel(mags_vec)/(time_win_data(i).Tend-time_win_data(i).Tstart);
HP(i).lamb=NaN; % numel(mags_vec(mags_vec>=mmin))/(time_win_data(i).Tend-time_win_data(i).Tstart);
switch iop
case 0
%OK
case 1
HP(i).lamb_all=HP(i).lamb_all*30;
HP(i).lamb=HP(i).lamb*30;
case 2
HP(i).lamb_all=HP(i).lamb_all*365;
HP(i).lamb=HP(i).lamb*365;
end
end
% K12NOV2015
end
cd ../
end
%%
function [MRPer,ExPr,MRPer_low,MRPer_high,ExPr_low,ExPr_high]=RetPeriodExcProbWrapper(lamb,CDF,DT,MaG)
DM=abs(CDF(:,1)-MaG);ind=find(DM==min(DM));
CDF1=CDF(ind,2);CDF_CI=[CDF(ind,3) CDF(ind,4)];
% Exceedance probability
EP=1-exp(-lamb*DT.*(1-CDF1));
EP_CI=1-exp(-lamb*DT.*(1-CDF_CI));
ExPr=EP;ExPr_low=EP_CI(:,2);ExPr_high=EP_CI(:,1);
% Mean Return Period
MRP=1/lamb./(1-CDF1);
MRP_CI=1/lamb./(1-CDF_CI);
MRPer=MRP;MRPer_low=MRP_CI(:,1);MRPer_high=MRP_CI(:,2);
end

99
SHAPE_Package/SHAPE_ver2b.0/SSH/ExcProbGRT.m Executable file
View File

@@ -0,0 +1,99 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [x,z]=ExcProbGRT(opt,xd,xu,dx,y,Mmin,lamb,eps,b,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dx<=0;error('Step must be greater than 0');end
%----------------------------------------------------------
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

78
SHAPE_Package/SHAPE_ver2b.0/SSH/ExcProbGRU.m Executable file
View File

@@ -0,0 +1,78 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [x,z]=ExcProbGRU(opt,xd,xu,dx,y,Mmin,lamb,eps,b)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dx<=0;error('Step must be greater than 0');end
%----------------------------------------------------------
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

116
SHAPE_Package/SHAPE_ver2b.0/SSH/ExcProbNPT.m Executable file
View File

@@ -0,0 +1,116 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [x,z]=...
ExcProbNPT(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dx<=0;error('Step must be greater than 0');end
%----------------------------------------------------------
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 %K15DEC2015
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

105
SHAPE_Package/SHAPE_ver2b.0/SSH/ExcProbNPU.m Executable file
View File

@@ -0,0 +1,105 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [x,z]=ExcProbNPU(opt,xd,xu,dx,y,Mmin,lamb,eps,h,xx,ambd)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dx<=0;error('Step must be greater than 0');end
%----------------------------------------------------------
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

59
SHAPE_Package/SHAPE_ver2b.0/SSH/Max_credM_GRT.m Executable file
View File

@@ -0,0 +1,59 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [T,m]=Max_credM_GRT(Td,Tu,dT,Mmin,lamb,eps,b,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dT<=0;error('Time Step must be greater than 0');end
%----------------------------------------------------------
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

63
SHAPE_Package/SHAPE_ver2b.0/SSH/Max_credM_GRU.m Executable file
View File

@@ -0,0 +1,63 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [T,m]=Max_credM_GRU(Td,Tu,dT,Mmin,lamb,eps,b)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dT<=0;error('Time Step must be greater than 0');end
%----------------------------------------------------------
T=(Td:dT:Tu)';
beta=b*log(10);
m=Mmin-eps/2+1/beta.*log(lamb*T);
end

98
SHAPE_Package/SHAPE_ver2b.0/SSH/Max_credM_NPT.m Executable file
View File

@@ -0,0 +1,98 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [T,m]=Max_credM_NPT(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dT<=0;error('Time Step must be greater than 0');end
%----------------------------------------------------------
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

98
SHAPE_Package/SHAPE_ver2b.0/SSH/Max_credM_NPT_O.m Executable file
View File

@@ -0,0 +1,98 @@
% [T,m]=Max_credM_NPT_O(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd,Mmax) ---- (Octave Compatible Version)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES THE MAXIMUM
% CREDIBLE MAGNITUDE VALUES FOR THE UPPER-BOUNDED NONPARAMETRIC
% DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [T,m]=Max_credM_NPT_O(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dT<=0;error('Time Step must be greater than 0');end
%----------------------------------------------------------
T=(Td:dT:Tu)';
n=length(T);
interval=[Mmin-eps/2 Mmax-0.001];
for i=1:n
m(i)=fzero(@(t) F_maxmagn(t,xx,h,ambd,Mmin-eps/2,Mmax,lamb,T(i)),interval);
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

98
SHAPE_Package/SHAPE_ver2b.0/SSH/Max_credM_NPU.m Executable file
View File

@@ -0,0 +1,98 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [T,m]=Max_credM_NPU(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dT<=0;error('Time Step must be greater than 0');end
%----------------------------------------------------------
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

99
SHAPE_Package/SHAPE_ver2b.0/SSH/Max_credM_NPU_O.m Executable file
View File

@@ -0,0 +1,99 @@
% [T,m]=Max_credM_NPU_O(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd) ---- (Octave Comlatible Version)
%
%USING THE NONPARAMETRIC ADAPTATIVE KERNEL APPROACH EVALUATES
% THE MAXIMUM CREDIBLE MAGNITUDE VALUES FOR THE UNBOUNDED
% NONPARAMETRIC DISTRIBUTION FOR MAGNITUDE.
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [T,m]=Max_credM_NPU_O(Td,Tu,dT,Mmin,lamb,eps,h,xx,ambd)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dT<=0;error('Time Step must be greater than 0');end
%----------------------------------------------------------
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));
m(i)=fzero(@(t) F_maxmagn_NPU(t,xx,h,ambd,Mmin-eps/2,lamb,T(i)),interval);
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

259
SHAPE_Package/SHAPE_ver2b.0/SSH/Nonpar.m Executable file
View File

@@ -0,0 +1,259 @@
% [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.
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UNBOUNDED
% NON-PARAMETRIC HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. The smoothing factor
% h, is estimated using the least-squares cross-validation for the Gaussian
% kernel function. The final form of the kernel is the adaptive kernel.
% In order to avoid repetitions, which cannot appear in a sample when the
% kernel estimators are used, the magnitude sample data are randomized
% within the magnitude round-off interval. The round-off interval length -
% eps is the least non-zero difference between sample data or 0.1 is the
% least difference if greater than 0.1. The randomization is done
% assuming exponential distribution of m in [m0-eps/2, m0+eps/2], where m0
% is the sample data point and eps is the length of roud-off inteval. The
% shape parameter of the exponential distribution is estimated from the whole
% data sample assuming the exponential distribution. The background sample
% - xx comprises the randomized values of magnitude doubled symmetrically
% with respect to the value Mmin-eps/2: length(xx)=2*length(M). Weigthing
% factors row vector for the adaptive kernel is of the same size as xx.
% 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.
%
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)); % K20OCT2014
lamb=30*nn/(t2-t1); % K20OCT2014
unit='month';
else
lamb_all=365*n/(t(n)-t(1)); % K20OCT2014
lamb=365*nn/(t2-t1); % K20OCT2014
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 = podwajanie(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 = podwajanie(x,x0)
% DOUBLES THE SAMPLE
% If the sample x(i) is is truncated from the left hand side and belongs
% to the interval [x0,inf) or it is truncated from the right hand side and
% belongs to the interval (-inf,x0]
% then the doubled sample is [-x(i)+2x0,x(i)]
% x - is the column data vector
% x2 - is the column vector of data doubled and sorted in the ascending
% order
x2=[-x+2*x0
x];
x2=sort(x2);
end
function [h,ierr]=hopt(x)
%Estimation of the optimal smoothing factor by means of the least squares
%method
% x - column data vector
% The result is an optimal smoothing factor
% ierr=0 - convergence, ierr=1 - multiple h, ierr=2 - approximate h is used
% The function calls the procedure FZERO for the function 'funct'
% NEW VERSION 2 - without a square matrix. Also equipped with extra zeros
% search
% MODIFIED JUNE 2014
ierr=0;
n=length(x);
x=sort(x);
interval=[0.000001 2*std(x)/n^0.2];
x1=funct(interval(1),x);
x2=funct(interval(2),x);
if x1*x2<0
[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

310
SHAPE_Package/SHAPE_ver2b.0/SSH/Nonpar_O.m Executable file
View File

@@ -0,0 +1,310 @@
% [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.
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UNBOUNDED
% NON-PARAMETRIC HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. The smoothing factor
% h, is estimated using the least-squares cross-validation for the Gaussian
% kernel function. The final form of the kernel is the adaptive kernel.
% In order to avoid repetitions, which cannot appear in a sample when the
% kernel estimators are used, the magnitude sample data are randomized
% within the magnitude round-off interval. The round-off interval length -
% eps is the least non-zero difference between sample data or 0.1 is the
% least difference if greater than 0.1. The randomization is done
% assuming exponential distribution of m in [m0-eps/2, m0+eps/2], where m0
% is the sample data point and eps is the length of roud-off inteval. The
% shape parameter of the exponential distribution is estimated from the whole
% data sample assuming the exponential distribution. The background sample
% - xx comprises the randomized values of magnitude doubled symmetrically
% with respect to the value Mmin-eps/2: length(xx)=2*length(M). Weigthing
% factors row vector for the adaptive kernel is of the same size as xx.
% 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.
%
function [lamb_all,lamb,lamb_err,unit,eps,ierr,h,xx,ambd]=...
Nonpar_O(t,M,iop,Mmin)
if isempty(t) || numel(t)<3 ,isempty(M(M>=Mmin)) %K03OCT
t=[1 2];M=[1 2]; end %K30SEP
lamb_err=0;
%%% %%%%%%%%%%%%%MICHAL
xx=NaN;
ambd=NaN;
%%% %%%%%%%%%%%%%MICHAL
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
% SL 03MAR2015 ----------------------------------
[NM,unit]=time_diff(t(1),t(n),iop);
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);
lamb=nn/NM;
% SL 03MAR2015 ----------------------------------
if nn<50
eps=0;ierr=0;h=0;
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 = podwajanie(xx,Mmin-eps/2);
[h,ierr]=hopt(xx);
[ambd]=scaling(xx,h);
% enai=dlmread('para.txt'); %for fixed xx,ambd to test in different platforms
% [ambd]=enai(:,1);
% xx=enai(:,2)';
% [h,ierr]=hopt(xx);
end
function [NM,unit]=time_diff(t1,t2,iop) % SL 03MAR2015
% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
%
% t1 - start time (in MATLAB numerical format)
% t2 - end time (in MATLAB numerical format) t2>=t1
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
% iop=2 - 'year'
%
% NM - number of time units from t1 to t2
% unit - string with name of time unit used ('year' or 'month' or 'day').
if iop==0
NM=(t2-t1);
unit='day';
elseif iop==1
V1=datevec(t1);
V2=datevec(t2);
NM=V2(3)/eomday(V2(1),V2(2))+V2(2)+12-V1(2)-V1(3)/eomday(V1(1),V1(2))...
+(V2(1)-V1(1)-1)*12;
unit='month';
else
V1=datevec(t1);
V2=datevec(t2);
NM2=V2(3);
if V2(2)>1
for k=1:V2(2)-1
NM2=NM2+eomday(V2(1),k);
end
end
day2=365; if eomday(V2(1),2)==29; day2=366; end;
NM2=NM2/day2;
NM1=V1(3);
if V1(2)>1
for k=1:V1(2)-1
NM1=NM1+eomday(V1(1),k);
end
end
day1=365; if eomday(V1(1),2)==29; day1=366; end;
NM1=(day1-NM1)/day1;
NM=NM2+NM1+V2(1)-V1(1)-1;
unit='year';
end
end
function [m_corr]=korekta(m,Mmin,eps,beta)
% RANDOMIZATION OF MAGNITUDE WITHIN THE ACCURACY INTERVAL
%
% m - input vector of magnitudes
% Mmin - catalog completeness level
% eps - accuracy of magnitude
% beta - the parameter of the unbounded exponential distribution
%
% m_corr - vector of randomized magnitudes
%
F1=1-exp(-beta*(m-Mmin-0.5*eps));
F2=1-exp(-beta*(m-Mmin+0.5*eps));
u=rand(size(m));
w=u.*(F2-F1)+F1;
m_corr=Mmin-log(1-w)./beta;
end
function x2 = podwajanie(x,x0)
% DOUBLES THE SAMPLE
% If the sample x(i) is is truncated from the left hand side and belongs
% to the interval [x0,inf) or it is truncated from the right hand side and
% belongs to the interval (-inf,x0]
% then the doubled sample is [-x(i)+2x0,x(i)]
% x - is the column data vector
% x2 - is the column vector of data doubled and sorted in the ascending
% order
x2=[-x+2*x0
x];
x2=sort(x2);
end
function [h,ierr]=hopt(x)
%Estimation of the optimal smoothing factor by means of the least squares
%method
% x - column data vector
% The result is an optimal smoothing factor
% ierr=0 - convergence, ierr=1 - multiple h, ierr=2 - approximate h is used
% The function calls the procedure FZERO for the function 'funct'
% NEW VERSION 2 - without a square matrix. Also equipped with extra zeros
% search
% MODIFIED JUNE 2014
ierr=0;
n=length(x);
x=sort(x);
interval=[0.000001 2*std(x)/n^0.2];
x1=funct(interval(1),x);
x2=funct(interval(2),x);
if x1*x2<0
fun = @(t) funct(t,x); % FOR OCTAVE
x0 =interval; % FOR OCTAVE
[hh(1),fval,exitflag] = fzero(fun,x0); % FOR OCTAVE
% Extra zeros search
jj=1;
for kk=2:7
interval(1)=1.1*hh(jj);
interval(2)=interval(1)+(kk-1)*hh(jj);
x1=funct(interval(1),x);
x2=funct(interval(2),x);
if x1*x2<0
jj=jj+1;
fun = @(t) funct(t,x); % FOR OCTAVE
x0 =interval; % FOR OCTAVE
[hh(jj),fval,exitflag] = fzero(fun,x0); % FOR OCTAVE
end
end
if jj>1;ierr=1;end
h=max(hh);
if exitflag==1;return;end
end
h=0.891836*(mean(x)-x(1))/(n^0.2);
ierr=2;
end
function [fct]=funct(t,x)
p2=1.41421356;
n=length(x);
yy=zeros(size(x));
for i=1:n,
xij=(x-x(i)).^2/t^2;
y=exp(-xij/4).*((xij/2-1)/p2)-2*exp(-xij/2).*(xij-1);
yy(i)=sum(y);
end;
fct=sum(yy)-2*n;
clear xij y yy;
end
function [ambd]=scaling(x,h)
% EVALUATES A VECTOR OF SCALING FACTORS FOR THE NONPARAMETRIC ADAPTATIVE
% ESTIMATION
% x - the n dimensional column vector of data values sorted in the ascending
% order
% h - the optimal smoothing factor
% ambd - the resultant n dimensional row vector of local scaling factors
n=length(x);
c=sqrt(2*pi);
gau=zeros(1,n);
for i=1:n,
gau(i)=sum(exp(-0.5*((x(i)-x)/h).^2))/c/n/h;
end
g=exp(mean(log(gau)));
ambd=sqrt(g./gau);
end
function [eps]=magn_accur(M)
x=sort(M);
d=x(2:length(x))-x(1:length(x)-1);
eps=min(d(d>0));
if eps>0.1; eps=0.1;end
end

373
SHAPE_Package/SHAPE_ver2b.0/SSH/Nonpar_tr.m Executable file
View File

@@ -0,0 +1,373 @@
% [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.
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UPPER-BOUNDED
% NON-PARAMETRIC HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% DESCRIPTION: The kernel estimator approach is a model-free alternative
% to estimating the magnitude distribution functions. The smoothing factor
% h, is estimated using the least-squares cross-validation for the Gaussian
% kernel function. The final form of the kernel is the adaptive kernel.
% In order to avoid repetitions, which cannot appear in a sample when the
% kernel estimators are used, the magnitude sample data are randomized
% within the magnitude round-off interval. The round-off interval length -
% eps is the least non-zero difference between sample data or 0.1 is the
% least difference if greater than 0.1. The randomization is done
% assuming exponential distribution of m in [m0-eps/2, m0+eps/2], where m0
% is the sample data point and eps is the length of roud-off inteval. The
% shape parameter of the exponential distribution is estimated from the whole
% data sample assuming the exponential distribution. The background sample
% - xx comprises the randomized values of magnitude doubled symmetrically
% with respect to the value Mmin-eps/2: length(xx)=2*length(M). Weigthing
% factors row vector for the adaptive kernel is of the same size as xx.
% The mean activity rate, lamb, is the number of events >=Mmin into the
% length of the period in which they occurred.
% The upper limit of the distribution Mmax is evaluated using
% the Kijko-Sellevol generic formula. If convergence is not reached the
% Whitlock @ Robson simplified formula is used:
% Mmaxest= 2(max obs M) - (second max obs M)).
%
% See: the references below for a more comprehensive description.
%
% This is a beta version of the program. Further developments are foreseen.
%
% REFERENCES:
%Silverman B.W. (1986) Density Estimation for Statistics and Data Analysis,
% Chapman and Hall, London
%Kijko A., Lasocki S., Graham G. (2001) Pure appl. geophys. 158, 1655-1665
%Lasocki S., Orlecka-Sikora B. (2008) Tectonophysics 456, 28-37
%Kijko, A., and M.A. Sellevoll (1989) Bull. Seismol. Soc. Am. 79, 3,645-654
%Lasocki, S., Urban, P. (2011) Acta Geophys 59, 659-673,
% doi: 10.2478/s11600-010-0049-y
%
% INPUT:
% t - vector of earthquake occurrence times
% M - vector of earthquake magnitudes (sample data)
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
% iop=2 - 'year'
% Mmin - lower bound of the distribution - catalog completeness level
%
% 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.
%
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)); % K20OCT2014
lamb=30*nn/(t2-t1); % K20OCT2014
unit='month';
else
lamb_all=365*n/(t(n)-t(1)); % K20OCT2014
lamb=365*nn/(t2-t1); % K20OCT2014
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 = podwajanie(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 = podwajanie(x,x0)
% DOUBLES THE SAMPLE
% If the sample x(i) is is truncated from the left hand side and belongs
% to the interval [x0,inf) or it is truncated from the right hand side and
% belongs to the interval (-inf,x0]
% then the doubled sample is [-x(i)+2x0,x(i)]
% x - is the column data vector
% x2 - is the column vector of data doubled and sorted in the ascending
% order
x2=[-x+2*x0
x];
x2=sort(x2);
end
function [h,ierr]=hopt(x)
%Estimation of the optimal smoothing factor by means of the least squares
%method
% x - column data vector
% The result is an optimal smoothing factor
% ierr=0 - convergence, ierr=1 - multiple h, ierr=2 - approximate h is used
% The function calls the procedure FZERO for the function 'funct'
% NEW VERSION 2 - without a square matrix. Also equipped with extra zeros
% search
% MODIFIED JUNE 2014
ierr=0;
n=length(x);
x=sort(x);
interval=[0.000001 2*std(x)/n^0.2];
x1=funct(interval(1),x);
x2=funct(interval(2),x);
if x1*x2<0
[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

506
SHAPE_Package/SHAPE_ver2b.0/SSH/Nonpar_tr_O.m Executable file
View File

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

83
SHAPE_Package/SHAPE_ver2b.0/SSH/Ret_periodGRT.m Executable file
View File

@@ -0,0 +1,83 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m,T]=Ret_periodGRT(Md,Mu,dM,Mmin,lamb,eps,b,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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

59
SHAPE_Package/SHAPE_ver2b.0/SSH/Ret_periodGRU.m Executable file
View File

@@ -0,0 +1,59 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m,T]=Ret_periodGRU(Md,Mu,dM,Mmin,lamb,eps,b)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
if Md<Mmin; Md=Mmin;end
m=(Md:dM:Mu)';
beta=b*log(10);
T=1/lamb./exp(-beta*(m-Mmin+eps/2));
end

94
SHAPE_Package/SHAPE_ver2b.0/SSH/Ret_periodNPT.m Executable file
View File

@@ -0,0 +1,94 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m,T]=Ret_periodNPT(Md,Mu,dM,Mmin,lamb,eps,h,xx,ambd,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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

91
SHAPE_Package/SHAPE_ver2b.0/SSH/Ret_periodNPU.m Executable file
View File

@@ -0,0 +1,91 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m,T]=Ret_periodNPU(Md,Mu,dM,Mmin,lamb,eps,h,xx,ambd)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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

250
SHAPE_Package/SHAPE_ver2b.0/SSH/TruncGR.m Executable file
View File

@@ -0,0 +1,250 @@
%
% [lamb_all,lamb,lmab_err,unit,eps,b,Mmax,err]=TruncGR(t,M,iop,Mmin)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% COMPLETE PART OF THE SAMPLE, 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
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UPPER-BOUNDED
% GUTENBERG-RICHETR HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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. Must be determined externally.
% 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.
%
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)); % K20OCT2014
lamb=30*nn/(t2-t1); % K20OCT2014
unit='month';
else
lamb_all=365*n/(t(n)-t(1)); % K20OCT2014
lamb=365*nn/(t2-t1); % K20OCT2014
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); %K21OCT2014
% x=sort(M,'descend');
% for i=1:n
% if x(i)<Mmin; break; end
% xx(i)=x(i); %
% end
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

387
SHAPE_Package/SHAPE_ver2b.0/SSH/TruncGR_O.m Executable file
View File

@@ -0,0 +1,387 @@
%
%[lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_Ob(t,M,iop,Mmin,Mmax,Nsynth)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE,
% THE GUTENBERG-RICHTER B-VALUE AND THE UPPER BOUND OF MAGNITUDE
% DISTRIBUTION USING THE UPPER-BOUNDED G-R LED MAGNITUDE DISTRIBUTION MODEL
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UPPER-BOUNDED
% GUTENBERG-RICHETR HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
%
% 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. Must be determined externally.
% Can take any value from [min(M), max(M)].
% Mmax - upper limit of Magnitude Distribution. Can be set by user, or
% estimate within the program - it then should be set as Mmax=[].
%
% OUTPUT:
%
% 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.
% BIAS - Mmax estimation Bias (Lasocki and Urban, 2011)
% SD - Mmax standard deviation (Lasocki ands Urban, 2011)
%
% LICENSE
% This file is a part of the IS-EPOS e-PLATFORM.
%
% This is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free
% Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
function [lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_O(t,M,iop,Mmin,Mmax,Nsynth)
if nargin==5;Nsynth=[];end % K08DEC2019
if isempty(t) || numel(t)<3 || isempty(M(M>=Mmin)) %K03OCT
t=[1 2];M=[1 2]; end %K30SEP
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
% SL 03MAR2015 ----------------------------------
[NM,unit]=time_diff(t(1),t(n),iop);
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);
lamb=nn/NM;
% SL 03MAR2015 ----------------------------------
if nn<15
eps=0;b=0;Mmax=0;err=0;
lamb_err=1;
BIAS=NaN;SD=NaN; %%% K 08NOV2019
return;
end
eps=magn_accur(M);
xx=M(M>=Mmin); %K21OCT2014
% x=sort(M,'descend');
% for i=1:n
% if x(i)<Mmin; break; end
% xx(i)=x(i); %
% end
clear x;
nn=length(xx);
Max_obs=max(xx);
beta0=0;
Mmax1=Max_obs;
if isempty(Mmax)==0 %%% K 28JUL2015
err_Mmax=0; %%% K 04DEC2019
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax); %%% K 28JUL2015
x0 = 1; %[0.05,4.0]; %%% K 28JUL2015 - See exception line 155
beta = fzero(fun,x0); %%% K 28JUL2015
err=0; %%% K 28JUL2015
else %%% K 28JUL2015 - line 150
err_Mmax=1; %%% K 04DEC2019
for i=1:50,
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax1);
x0 =1; %[0.05,4.0]; %%% K29JUL2015 - See exception line 155
beta = fzero(fun,x0);
Mmax=Max_obs+moja_calka('f_podc',Mmin,Max_obs,1e-5,nn,beta,Mmin-eps/2,Mmax1);
if ((abs(Mmax-Mmax1)<0.01)&&(abs(beta-beta0)<0.0001))
err=0;
break;
end
Mmax1=Mmax;
beta0=beta;
end
if i==50;
err=1.0;
Mmax=2*xx(1)-xx(2);
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax);
x0 =1;
beta = fzero(fun,x0);
end
end %%% K 28JUL2015
b=beta/log(10);
clear xx
% Exception for v-value
if b<0.05 || b>6.0; error('Unacceptable b-value, abort and select different dataset');end
beta;
% Estimation of Mmax Bias %%% K 04DEC2019
% (Lasocki and Urban, 2011, doi:10.2478/s11600-010-0049-y)
if isempty(Nsynth)==0 && err_Mmax==1 % set number of synthetic datasets, e.g. 10000
[BIAS,SD]=Mmax_Bias_GR(t,M,Mmin,Mmax,b,err,Nsynth);
elseif isempty(Nsynth)==0 && err_Mmax==0;
warning('Mmax must be empty for BIAS calculation');BIAS=[];SD=[];
else; BIAS=0;SD=0;
end
end
function [NM,unit]=time_diff(t1,t2,iop) % SL 03MAR2015
% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
%
% 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 [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
% --------------------- Mmax BIAS estimation routine ---------------------- K 08NOV2019
function [BIAS,SD]=Mmax_Bias_GR(t,m,Mc,Mmax1,b,err,synth)
if err~=0
warning('process did not converge!!'),pause
end
MAXm=max(m);beta=b*log(10);N=numel(m(m>=Mc));DeltaM=MAXm-Mc;
for j=1:synth %set number of synthetic datasets, default is 10000
% % CDF:
M=Mc:0.0001:MAXm;upt=1-exp(-beta*(M-Mc));
dwt=1-exp(-beta*(MAXm-Mc));F=upt./dwt; % j
% linear interpolation to assign magnitude values from a uniform distribution sample
iM=rand(1,N);M1=interp1q(F',M',iM');
br(j)=1/(log(10)*(mean(M1)-min(M1)));DM=range(M1);
Mmax=max(M1);
% Iteration Process to estimate b and Mmax
b1=1;best=[1.0 10.0];i=1;
while min(abs(diff(best)))>0.00001
w=exp(b1*(Mmax-Mc));E1=expint(N/(w-1));E2=expint(N*w/(w-1));
%E=expint(w);
Mme=Mmax+(E1-E2)/(b1*exp(-N/(w-1)))+(Mc)*exp(-N); %Mme=round(Mme/EPS)*EPS;
if isnan(Mme)
KM=sort(unique(M1),'descend');
Mme=2*KM(1)-KM(2);
end
fun=@(bb) 1/bb+(Mme-Mc)/(1-exp(bb*(Mme-Mc)))-mean(M1)+Mc; %consider th5 last 0.05 term
b1=fzero(fun,1);best(i)=b1;i=i+1;
if i==50
warning('process did not converge!!');break
end
end
be(j)=b1/log(10);
Me(j)=Mme;dm(j)=DM;Mm(j)=Mmax;
end
BIAS=mean(MAXm-Me)
SD=std(MAXm-Me);
%b-mean(be) %check b-value difference
%histogram(be)
% MAXm: maximum magnitude in the real catalog
% Mmax: maximum magnitudes observed in the synthetic catalogs (rounded)
% Me: maximum magnitude estimates for the synthetic catalogs
% Mmax1: maximum magnitude estimated by GRT
end

387
SHAPE_Package/SHAPE_ver2b.0/SSH/TruncGR_O_CI.m Executable file
View File

@@ -0,0 +1,387 @@
%
%[lamb_all,lamb,lamb_err,unit,eps,b,Mmax,err,BIAS,SD]=TruncGR_Ob(t,M,iop,Mmin,Mmax,Nsynth)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE,
% THE GUTENBERG-RICHTER B-VALUE AND THE UPPER BOUND OF MAGNITUDE
% DISTRIBUTION USING THE UPPER-BOUNDED G-R LED MAGNITUDE DISTRIBUTION MODEL
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UPPER-BOUNDED
% GUTENBERG-RICHETR HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
%
% 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. Must be determined externally.
% Can take any value from [min(M), max(M)].
% Mmax - upper limit of Magnitude Distribution. Can be set by user, or
% estimate within the program - it then should be set as Mmax=[].
%
% OUTPUT:
%
% 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.
% BIAS - Mmax estimation Bias (Lasocki and Urban, 2011)
% SD - Mmax standard deviation (Lasocki ands Urban, 2011)
%
% LICENSE
% This file is a part of the IS-EPOS e-PLATFORM.
%
% This is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free
% Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
function [b,lamb_all,lamb,lamb_err,unit,eps,Mmax,err,BIAS,SD]=TruncGR_O_CI(t,M,iop,Mmin,Mmax,Nsynth)
if nargin==5;Nsynth=[];end % K08DEC2019
if isempty(t) || numel(t)<3 || isempty(M(M>=Mmin)) %K03OCT
t=[1 2];M=[1 2]; end %K30SEP
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
% SL 03MAR2015 ----------------------------------
[NM,unit]=time_diff(t(1),t(n),iop);
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);
lamb=nn/NM;
% SL 03MAR2015 ----------------------------------
if nn<15
eps=0;b=0;Mmax=0;err=0;
lamb_err=1;
BIAS=NaN;SD=NaN; %%% K 08NOV2019
return;
end
eps=magn_accur(M);
xx=M(M>=Mmin); %K21OCT2014
% x=sort(M,'descend');
% for i=1:n
% if x(i)<Mmin; break; end
% xx(i)=x(i); %
% end
clear x;
nn=length(xx);
Max_obs=max(xx);
beta0=0;
Mmax1=Max_obs;
if isempty(Mmax)==0 %%% K 28JUL2015
err_Mmax=0; %%% K 04DEC2019
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax); %%% K 28JUL2015
x0 = 1; %[0.05,4.0]; %%% K 28JUL2015 - See exception line 155
beta = fzero(fun,x0); %%% K 28JUL2015
err=0; %%% K 28JUL2015
else %%% K 28JUL2015 - line 150
err_Mmax=1; %%% K 04DEC2019
for i=1:50,
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax1);
x0 =1; %[0.05,4.0]; %%% K29JUL2015 - See exception line 155
beta = fzero(fun,x0);
Mmax=Max_obs+moja_calka('f_podc',Mmin,Max_obs,1e-5,nn,beta,Mmin-eps/2,Mmax1);
if ((abs(Mmax-Mmax1)<0.01)&&(abs(beta-beta0)<0.0001))
err=0;
break;
end
Mmax1=Mmax;
beta0=beta;
end
if i==50;
err=1.0;
Mmax=2*xx(1)-xx(2);
fun = @(b) bet_est(b,mean(xx),Mmin-eps/2,Mmax);
x0 =1;
beta = fzero(fun,x0);
end
end %%% K 28JUL2015
b=beta/log(10);
clear xx
% Exception for v-value
if b<0.05 || b>6.0; error('Unacceptable b-value, abort and select different dataset');end
beta;
% Estimation of Mmax Bias %%% K 04DEC2019
% (Lasocki and Urban, 2011, doi:10.2478/s11600-010-0049-y)
if isempty(Nsynth)==0 && err_Mmax==1 % set number of synthetic datasets, e.g. 10000
[BIAS,SD]=Mmax_Bias_GR(t,M,Mmin,Mmax,b,err,Nsynth);
elseif isempty(Nsynth)==0 && err_Mmax==0;
warning('Mmax must be empty for BIAS calculation');BIAS=[];SD=[];
else; BIAS=0;SD=0;
end
end
function [NM,unit]=time_diff(t1,t2,iop) % SL 03MAR2015
% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
%
% 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 [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
% --------------------- Mmax BIAS estimation routine ---------------------- K 08NOV2019
function [BIAS,SD]=Mmax_Bias_GR(t,m,Mc,Mmax1,b,err,synth)
if err~=0
warning('process did not converge!!'),pause
end
MAXm=max(m);beta=b*log(10);N=numel(m(m>=Mc));DeltaM=MAXm-Mc;
for j=1:synth %set number of synthetic datasets, default is 10000
% % CDF:
M=Mc:0.0001:MAXm;upt=1-exp(-beta*(M-Mc));
dwt=1-exp(-beta*(MAXm-Mc));F=upt./dwt; % j
% linear interpolation to assign magnitude values from a uniform distribution sample
iM=rand(1,N);M1=interp1q(F',M',iM');
br(j)=1/(log(10)*(mean(M1)-min(M1)));DM=range(M1);
Mmax=max(M1);
% Iteration Process to estimate b and Mmax
b1=1;best=[1.0 10.0];i=1;
while min(abs(diff(best)))>0.00001
w=exp(b1*(Mmax-Mc));E1=expint(N/(w-1));E2=expint(N*w/(w-1));
%E=expint(w);
Mme=Mmax+(E1-E2)/(b1*exp(-N/(w-1)))+(Mc)*exp(-N); %Mme=round(Mme/EPS)*EPS;
if isnan(Mme)
KM=sort(unique(M1),'descend');
Mme=2*KM(1)-KM(2);
end
fun=@(bb) 1/bb+(Mme-Mc)/(1-exp(bb*(Mme-Mc)))-mean(M1)+Mc; %consider th5 last 0.05 term
b1=fzero(fun,1);best(i)=b1;i=i+1;
if i==50
warning('process did not converge!!');break
end
end
be(j)=b1/log(10);
Me(j)=Mme;dm(j)=DM;Mm(j)=Mmax;
end
BIAS=mean(MAXm-Me)
SD=std(MAXm-Me);
%b-mean(be) %check b-value difference
%histogram(be)
% MAXm: maximum magnitude in the real catalog
% Mmax: maximum magnitudes observed in the synthetic catalogs (rounded)
% Me: maximum magnitude estimates for the synthetic catalogs
% Mmax1: maximum magnitude estimated by GRT
end

167
SHAPE_Package/SHAPE_ver2b.0/SSH/UnlimitGR.m Executable file
View File

@@ -0,0 +1,167 @@
% [lamb_all,lamb,lmab_err,unit,eps,b]=UnlimitGR(t,M,iop,Mmin)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE AND THE
% GUTENBERG-RICHTER B-VALUE USING THE UNLIMITED G-R LED MAGNITUDE
% DISTRIBUTION MODEL
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UNBOUNDED
% GUTENBERG-RICHETR HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
%
% 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 if the least non-zero difference
% between sample data or 0.1 is 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. Must be determined externally.
% 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.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
function [lamb_all,lamb,lamb_err,unit,eps,b,bCI]=UnlimitGR(t,M,iop,Mmin,Nbst)
if isempty(t) || numel(t)<3 || isempty(M(M>=Mmin)) %K03OCT
t=[1 2];M=[1 2]; end %K30SEP
lamb_err=0;
n=length(M);
t1=t(1);
for i=1:n
if M(i)>=Mmin; break; end
t1=t(i+1);
end
t2=t(n);
for i=n:1
if M(i)>=Mmin; break; end
t2=t(i-1);
end
nn=0;
for i=1:n
if M(i)>=Mmin
nn=nn+1;
end
end
% SL 03MAR2015 ----------------------------------
[NM,unit]=time_diff(t(1),t(n),iop);
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);
lamb=nn/NM;
% SL 03MAR2015 ----------------------------------
if nn<7
eps=0;b=0;
lamb_err=1;
return;
end
eps=magn_accur(M);
xx=M(M>=Mmin); %K21OCT2014
% x=sort(M,'descend');
% for i=1:n
% if x(i)<Mmin; break; end
% xx(i)=x(i); %
% end
clear x;
beta=1/(mean(xx)-Mmin+eps/2);
b=beta/log(10);
% Add bootstrap CI of b-value %K03JUN2020
b_aki=@(xx)1/(log(10)*(mean(xx)-min(xx)+eps/2)); %K03JUN2020
bCI=bootci(Nbst,{b_aki,xx},'alpha',0.05); %K03JUN2020
clear xx
end
function [NM,unit]=time_diff(t1,t2,iop) % SL 03MAR2015
% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
%
% t1 - start time (in MATLAB numerical format)
% t2 - end time (in MATLAB numerical format) t2>=t1
% iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month',
% iop=2 - 'year'
%
% NM - number of time units from t1 to t2
% unit - string with name of time unit used ('year' or 'month' or 'day').
if iop==0
NM=(t2-t1);
unit='day';
elseif iop==1
V1=datevec(t1);
V2=datevec(t2);
NM=V2(3)/eomday(V2(1),V2(2))+V2(2)+12-V1(2)-V1(3)/eomday(V1(1),V1(2))...
+(V2(1)-V1(1)-1)*12;
unit='month';
else
V1=datevec(t1);
V2=datevec(t2);
NM2=V2(3);
if V2(2)>1
for k=1:V2(2)-1
NM2=NM2+eomday(V2(1),k);
end
end
day2=365; if eomday(V2(1),2)==29; day2=366; end;
NM2=NM2/day2;
NM1=V1(3);
if V1(2)>1
for k=1:V1(2)-1
NM1=NM1+eomday(V1(1),k);
end
end
day1=365; if eomday(V1(1),2)==29; day1=366; end;
NM1=(day1-NM1)/day1;
NM=NM2+NM1+V2(1)-V1(1)-1;
unit='year';
end
end
function [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

64
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_GRT.m Executable file
View File

@@ -0,0 +1,64 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m, PDF_GRT, CDF_GRT]=dist_GRT(Md,Mu,dM,Mmin,eps,b,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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

69
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_GRT_CI.m Executable file
View File

@@ -0,0 +1,69 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m, PDF_GRT, CDF_GRT]=dist_GRT_CI(Md,Mu,dM,Mmin,eps,b,Mmax,bCI)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
m=(Md:dM:Mu)';
beta=b*log(10);beta_low=bCI(1)*log(10);beta_high=bCI(2)*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));
CDF_low=1-exp(-beta_low*(m-Mmin+eps/2));
CDF_high=1-exp(-beta_high*(m-Mmin+eps/2));
CDF_GRT=[CDF_GRT CDF_low CDF_high];
end

61
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_GRU.m Executable file
View File

@@ -0,0 +1,61 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m, PDF_GRU, CDF_GRU]=dist_GRU(Md,Mu,dM,Mmin,eps,b)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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

67
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_GRU_CI.m Executable file
View File

@@ -0,0 +1,67 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m, PDF_GRU, CDF_GRU]=dist_GRU_CI(Md,Mu,dM,Mmin,eps,b,bCI)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
m=(Md:dM:Mu)';
beta=b*log(10);beta_low=bCI(1)*log(10);beta_high=bCI(2)*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));
CDF_low=1-exp(-beta_low*(m-Mmin+eps/2));
CDF_high=1-exp(-beta_high*(m-Mmin+eps/2));
CDF_GRU=[CDF_GRU CDF_low CDF_high];
end

116
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_NPT.m Executable file
View File

@@ -0,0 +1,116 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m,PDF_NPT,CDF_NPT]=dist_NPT(Md,Mu,dM,Mmin,eps,h,xx,ambd,Mmax)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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. Use
% "podwajanie.m" first to accmoplish that.
% 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

116
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_NPT_CI.m Executable file
View File

@@ -0,0 +1,116 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [CDF_NPT]=dist_NPT_CI(Md,Mu,dM,Mmin,eps,h,data,Mmax)
xx=data(:,1);ambd=data(:,2)';
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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. Use
% "podwajanie.m" first to accmoplish that.
% 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

114
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_NPU.m Executable file
View File

@@ -0,0 +1,114 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [m, PDF_NPU, CDF_NPU]=dist_NPU(Md,Mu,dM,Mmin,eps,h,xx,ambd)
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
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. Use
% "podwajanie.m" first to accmoplish that.
% 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

118
SHAPE_Package/SHAPE_ver2b.0/SSH/dist_NPU_CI.m Executable file
View File

@@ -0,0 +1,118 @@
% [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: S. Lasocki 06/2014 within IS-EPOS project.
%
% 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.
%
function [CDF_NPU]=dist_NPU_CI(Md,Mu,dM,Mmin,eps,h,data)
xx=data(:,1);ambd=data(:,2)';
% -------------- VALIDATION RULES ------------- K_21NOV2016
if dM<=0;error('Magnitude Step must be greater than 0');end
%----------------------------------------------------------
m=(Md:dM:Mu)';
nn=length(m);
Mmax=10000;
mian=2*(Dystr_npr(Mmax,xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h));
% mian=1;
for i=1:nn
if m(i)>=Mmin-eps/2
PDF_NPU(i)=dens_npr1(m(i),xx,ambd,h,Mmin-eps/2)/mian;
CDF_NPU(i)=2*(Dystr_npr(m(i),xx,ambd,h)-Dystr_npr(Mmin-eps/2,xx,ambd,h))/mian;
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. Use
% "podwajanie.m" first to accmoplish that.
% 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

5
SHAPE_Package/SHAPE_ver2b.0/TIME_WINDOWS/ST2_test_timewindows.txt Executable file
View File

@@ -0,0 +1,5 @@
735470.75 735679.15
735679.15 735808.95
735808.95 735968.08
735968.08 736027.27
736027.27 736191.58

148
SHAPE_Package/SHAPE_ver2b.0/Zplo_ver2b_1.m Executable file
View File

@@ -0,0 +1,148 @@
close all;d=figure('Position',[25 00 1500 800]);
% check whether the selected time windows are overlapping or not
TT=[];Tcat=Catalog(1).val;Ncat=Tcat(Tcat>=time_windows(1).Tstart & Tcat<=time_windows(length(ExPr)).Tend);
for i=1:length(MRPer);
TW1(i)=time_windows(i).Tstart;Tw2(i)=time_windows(i).Tend;
tplo(i)=mean([time_windows(i).Tstart time_windows(i).Tend]);meanM(i)=mean(time_windows(i).M);hold on
TT=[TT;time_windows(i).Time];
lambda(i)=HP(i).lamb;
if strcmp(HP(1).method,'GRU') || strcmp(HP(1).method,'GRT');yyaxis right;bval(i)=HP(i).b;end
end
if (strcmp(Plotopt,'ON'))
%if numel(TT)==numel(Ncat)
DTW=TW1(2:length(TW1))-Tw2(1:length(Tw2)-1); %%% THIS SEEMS TO WORK!!!!
if isempty(find(DTW<0))
overlap='NO';
for i=1:length(MRPer);
subplot(3,1,1) % plot Mean return period
hold on;fill([time_windows(i).Tstart time_windows(i).Tend time_windows(i).Tend time_windows(i).Tstart],...
[MRPer_high(i) MRPer_high(i) MRPer_low(i) MRPer_low(i)],[0.91 0.91 0.91],'facealpha',0.25)
plot([time_windows(i).Tstart time_windows(i).Tend],[MRPer(i) MRPer(i)],'k-','LineWidth',2)
if i<length(MRPer);plot([time_windows(i).Tend time_windows(i+1).Tstart],[MRPer(i) MRPer(i+1)],'k--');end
datetick('x',20);title(['Mean Return Period for M\geq',num2str(MaG)],'FontSize',16);ylabel([Tunit,'s'],'FontSize',18)
subplot(3,1,2) % plot Exceedance Probability
hold on;fill([time_windows(i).Tstart time_windows(i).Tend time_windows(i).Tend time_windows(i).Tstart],...
[ExPr_high(i) ExPr_high(i) ExPr_low(i) ExPr_low(i)],[0.91 0.91 0.91],'facealpha',0.25)
plot([time_windows(i).Tstart time_windows(i).Tend],[ExPr(i) ExPr(i)],'k-','LineWidth',2)
if i<length(ExPr);plot([time_windows(i).Tend time_windows(i+1).Tstart],[ExPr(i) ExPr(i+1)],'k--');end
datetick('x',20);title(['Exceedance Probability for M\geq',num2str(MaG),' within ',num2str(Plength), ' ',Tunit,'(s) period'],'FontSize',16);ylabel('probability','FontSize',14)
subplot(3,1,3) % plot Activity rate
hold on;yyaxis left;plot([time_windows(i).Tstart time_windows(i).Tend],[HP(i).lamb HP(i).lamb],'k-','LineWidth',2)
if i<length(ExPr);plot([time_windows(i).Tend time_windows(i+1).Tstart],[HP(i).lamb HP(i+1).lamb],'k--');end
datetick('x',20);title(['Activity Rate'],'FontSize',16);ylabel(['Events/',Tunit],'FontSize',14,'Color','k')
set(gca,'YColor','k');
% plot b-value (GR) or mean M (NP)
if strcmp(HP(1).method,'GRU') || strcmp(HP(1).method,'GRT');yyaxis right;
fill([time_windows(i).Tstart time_windows(i).Tend time_windows(i).Tend time_windows(i).Tstart],...
[HP(i).bCI(2) HP(i).bCI(2) HP(i).bCI(1) HP(i).bCI(1)],[0.99 0.81 0.31],'LineStyle','-','Marker','none','facealpha',0.25)
plot([time_windows(i).Tstart time_windows(i).Tend],[HP(i).b HP(i).b],'-','LineWidth',2)
ylabel('b-value','FontSize',14);
else
yyaxis right;plot([time_windows(i).Tstart time_windows(i).Tend],[mean(time_windows(i).M) mean(time_windows(i).M)],'-','LineWidth',2)
ylabel('mean Magnitude','FontSize',14);
end
end
else
overlap='YES';
subplot(3,1,1) % plot Mean return period
errorbar(tplo,MRPer,MRPer-MRPer_low,MRPer_high-MRPer,'o','LineWidth',1,'MarkerSize',8);hold on;
plot(tplo,MRPer,'ko','LineWidth',2,'MarkerSize',4);
datetick('x',20);title(['Mean Return Period for M\geq',num2str(MaG)],'FontSize',16);ylabel([Tunit,'s'],'FontSize',18)
subplot(3,1,2) % plot Exceedance Probability
errorbar(tplo,ExPr,ExPr-ExPr_low,ExPr_high-ExPr,'o','LineWidth',1,'MarkerSize',8);hold on
plot(tplo,ExPr,'ko','LineWidth',2,'MarkerSize',4);
datetick('x',20);title(['Exceedance Probability for M\geq',num2str(MaG),' within ',num2str(Plength), ' ',Tunit,'(s) period'],'FontSize',16);ylabel('probability','FontSize',14)
subplot(3,1,3) % plot Activity rate
plot(tplo,lambda,'o','LineWidth',2,'MarkerSize',8);ylabel(['Events/',Tunit],'FontSize',14)
if strcmp(HP(1).method,'GRU') || strcmp(HP(1).method,'GRT');
yyaxis right;
for i=1:length(HP);bval_low(i)=HP(i).bCI(1);bval_high(i)=HP(i).bCI(2);end
errorbar(tplo,bval,bval-bval_low,bval_high-bval,'o','LineWidth',1,'MarkerSize',8);hold on
plot(tplo,bval,'ko','LineWidth',2,'MarkerSize',4);
ylabel('b-value','FontSize',14);
else; yyaxis right;plot(tplo,meanM,'o','LineWidth',2,'MarkerSize',8)
ylabel('mean Magnitude','FontSize',14);end
datetick('x',20);title('Activity Rate','FontSize',16);
end
if isempty(PROD_Data)==0
subplot(3,1,1);yyaxis right;plot(PROD_Data(1).val,PROD_Data(PROD_FIELD).val,'-','Linewidth',1);ylabel(PROD_Data(PROD_FIELD).field,'interpreter','none','FontSize',14);
subplot(3,1,2);yyaxis right;plot(PROD_Data(1).val,PROD_Data(PROD_FIELD).val,'-','Linewidth',1);ylabel(PROD_Data(PROD_FIELD).field,'interpreter','none','FontSize',14);
end
subplot(3,1,3);xlabel('Date','FontSize',18)
% option to switch linear-log Y axis Scale
txt = uicontrol('Parent',d,...
'Style','text',...
'Position',[200 470 150 30],...
'String','Select Y Axis Scale:');
popup = uicontrol('Parent',d,...
'Style','popup',...
'Position',[350 480 120 25],...
'String',{'Linear';'Log'},...
'Callback',@popup_callback);
btn = uicontrol('Parent',d,...
'Position',[210 688 210 50],...
'String','SAVE and CLOSE',...
'FontSize',18,...
'ForeGroundColor','r',...
'FontWeight','Bold',...
'Callback',@savefig_callback);
choice = 'Linear';
% Wait for d to close before running to completion
uiwait(d);
elseif (strcmp(Plotopt,'OFF'));close all
if numel(TT)==numel(Ncat)
overlap='NO';else; overlap='YES';end
end
function popup_callback(popup,event)
idx = popup.Value;
popup_items = popup.String;
% This code uses dot notation to get properties.
% Dot notation runs in R2014b and later.
% For R2014a and earlier:
% idx = get(popup,'Value');
% popup_items = get(popup,'String');
choice = char(popup_items(idx,:));
subplot(3,1,1);yyaxis left;
set(gca,'YScale',choice)
end
function savefig_callback(popup,event)
cd Outputs_SHA
print(gcf,'SHA.jpeg','-djpeg','-r300')
savefig(gcf,'SHA.fig')
% This code uses dot notation to get properties.
% Dot notation runs in R2014b and later.
% For R2014a and earlier:
% idx = get(popup,'Value');
% popup_items = get(popup,'String');
cd ../
delete(gcf)
end

93
SHAPE_Package/SHAPE_ver2b.0/Zsave_output_ver2b_1.m Executable file
View File

@@ -0,0 +1,93 @@
% ---- Save *.txt file with Parameters Report ----
cd Outputs_SHA
fid=fopen('REPORT_Hazard_Analysis.txt','w');
fprintf(fid,['Parameters Report & Results for HAZARD ANALYSIS (created on ', datestr(now),')\n']);
fprintf(fid,['Parameters Estimated: Mean Return Period (MRP) and Exceedance Probability (EP) \n']);
fprintf(fid,'------------------------------------------------------------------------------------\n');
fprintf(fid,['<Magnitude Scale Selected >: ', MScale,'\n']);
fprintf(fid,['<Time Unit >: ', Tunit,'\n']);
fprintf(fid,['<Magnitude Range >: ', num2str(Mc), ' to ', num2str(max(Cmag)),'\n']);
if strcmp(method,'GRT')==1 || strcmp(method,'NPT')==1
fprintf(fid,['<Maximum Magnitude >: ', num2str(HP(1).Mmax),'\n']);
else
fprintf(fid,['<Maximum Magnitude >: Unbounded','\n']);
end
fprintf(fid,['<Magnitude Distribution Model >: ', method,'\n']);
fprintf(fid,['<Magnitude (for EPP and MRP) >: ', num2str(MaG),'\n']);
fprintf(fid,['<Time Period (for EPR) >: ', num2str(Plength),' ',Tunit,'s','\n']);
fprintf(fid,['<Time Window Creation Mode >: ', winmode,'\n']);
if strcmp(winmode,'Time')==1
fprintf(fid,['< Window Size >: ', num2str(window_size),'(days) \n']);
fprintf(fid,['< Window Step >: ', num2str(dt),'(days) \n']);
elseif strcmp(winmode,'Events')==1
fprintf(fid,['< Window Size >: ', num2str(window_size),'(events) \n']);
fprintf(fid,['< Window Step >: ', num2str(dt),'(days) \n']);
elseif strcmp(winmode,'Graphical')==1
fprintf(fid,['< Window Size >: variable \n']);
fprintf(fid,['< Window Step >: variable \n']);
end
fprintf(fid,['<Overlapping Time Windows >: ', overlap,'\n']);
fprintf(fid,['<Bootstrap Iterations >: ', num2str(Nbst),'\n']);
fprintf(fid,'------------------------------------------------------------------------------------\n');
for j=1:numel(HP)
SN(j)=j;Nevents(j)=numel(time_windows(j).M);TS(j)=time_windows(j).Tstart;TE(j)=time_windows(j).Tend;
end
fprintf(fid,[' Set N Starting Date/Time Ending Date/Time events MRP MRP0.95CI EP E)0.95CI b-value b-value0.95CI \n']);
fprintf(fid,[' per ',Tunit, ' ',Tunit,'s' '\n']);
for i=1:numel(HP)
if strcmp(method,'GRU')==1 || strcmp(method,'GRT')==1;
fprintf(fid,['%4d %5d %s %s %9.3f %13.3f %13.3f - %13.3f %13.11f %13.11f - %13.11f %5.3f %5.3f %5.3f \n'],SN(i),Nevents(i),datestr(TS(i),0),datestr(TE(i),0),lambda(i),MRPer(i),...
MRPer_low(i),MRPer_high(i),ExPr(i),ExPr_low(i),ExPr_high(i),bval(i),HP(i).bCI(1),HP(i).bCI(2));
else
fprintf(fid,['%4d %5d %s %s %9.3f %13.3f %13.3f - %13.3f %13.11f %13.3f - %13.11f %s %s \n'],SN(i),Nevents(i),datestr(TS(i),0),datestr(TE(i),0),lambda(i),MRPer(i),...
MRPer_low(i),MRPer_high(i),ExPr(i),ExPr_low(i),ExPr_high(i),'NaN','NaN');
end
end
fclose(fid);
% Save output structure time_window merged with HP
for i=1:length(HP)
SHA(i).Time=time_windows(i).Time;
SHA(i).M=time_windows(i).M;
SHA(i).Mmin=HP(i).mmin;
SHA(i).eps=HP(i).eps;
SHA(i).lambd=HP(i).lamb;
SHA(i).lambd_err=HP(i).lamb_err;
SHA(i).unit=HP(i).unit;
SHA(i).method=HP(i).method;
if strcmp(method,'GRU')==1 || strcmp(method,'GRT')==1
SHA(i).b=HP(i).b;
SHA(i).b_0025=HP(i).bCI(1);
SHA(i).b_0975=HP(i).bCI(2);
else
SHA(i).h=HP(i).h;
SHA(i).xx=HP(i).xx;
SHA(i).ambd=HP(i).ambd;
SHA(i).ierr=HP(i).ierr;
end
if strcmp(method,'GRT')==1 || strcmp(method,'NPT')==1
SHA(i).Mmax=HP(i).Mmax;
SHA(i).err=HP(i).err;
else
end
SHA(i).PDF=HP(i).PDF; %K29JAN2020
SHA(i).CDF=HP(i).CDF; %K29JAN2020
SHA(i).MRP=MRPer(i); %K03JUN2020
SHA(i).MRP_0025=MRPer_low(i); %K03JUN2020
SHA(i).MRP_0975=MRPer_high(i); %K03JUN2020
SHA(i).EP_=ExPr(i); %K03JUN2020
SHA(i).EP_0025=ExPr_low(i); %K03JUN2020
SHA(i).EP_0975=ExPr_high(i); %K03JUN2020
end
save('SHA.mat','SHA')
cd ../