magnitude_det_rate_estimation_mle.ipynb

_magnitude_det_rate_estimation_mle.ipynb

catalog.txt

1966 09 14 05 42 29	51.6	143.9	15	3.6
1967 01 09 20 26 09	51.5	143.5	10	3.6
1967 05 21 03 52 25	51.4	144	15	3.4
1968 01 25 23 53 50	51.3	143.4	15	3.6
1970 11 04 11 34 53	51.2	143.4	15	3.3
1970 11 14 17 48 45	51.4	143.1	15	3.4
1972 09 09 23 06 11	51.2	143.5	10	3.3
1975 01 04 03 58 48	51.4	143.5	15	3.9
1976 04 16 08 30 33	51.1	143.3	10	3.9
1977 10 23 05 07 35	51.6	143	15	3.4
1988 07 23 19 31 20	51.46	143.81	7	3.5
1988 08 24 05 44 49	51.36	143.48	9	4.5
1989 04 17 03 50 34	51.55	143.01	8	4.2
1992 01 06 14 21 33	51.14	143.36	9	3.1
1994 04 06 05 40 16	51.3	143.8	6	3.1
1995 10 05 19 27 26	51.35	143.45	10	3
1997 01 22 20 51 23	51.32	143.61	10	3.6
1997 07 09 08 00 43	51.56	143.28	9	3
1997 07 10 22 17 02	51.54	143.03	10	3.3
1998 04 30 21 31 39	51.43	143	9	3
2000 05 09 05 03 16	51.59	144	9	3.2
2002 01 03 10 49 58	51.22	143.83	14	3.6
2006 10 03 14 37 19	51.422	143.058	4.1	1.4
2006 10 13 20 04 11	51.569	143.085	10	3.2
2006 11 01 19 36 08	51.444	143.2	13.7	1.2
2006 11 28 07 09 50	51.106	143.014	18.3	1.3
2007 02 15 12 15 23	51.327	143.369	4.1	1.9
2007 02 20 17 14 18	51.28	143.55	2	1.5
2007 03 16 08 41 52	51.494	143.008	0	0.7
2007 04 01 17 57 10	51.531	143.16	0	2.1
2007 04 05 16 41 29	51.458	143.114	24.2	1.4
2007 06 12 10 49 58	51.397	143.039	0	1.2
2007 06 16 13 44 23	51.272	143.374	26.5	1.1
2007 06 25 15 57 01	51.469	143.111	5.2	1.9
2007 08 25 04 51 10	51.287	143.368	23.7	2.5
2007 08 31 12 02 32	51.237	143.339	3	1.4
2007 09 19 17 34 21	51.357	143.084	4.2	1.4
2008 08 08 20 35 25	51.29	143.035	10	0.6
2009 02 24 10 21 12	51.229	143.166	12.1	1.1
2009 04 01 12 29 50	51.175	143.586	2	1.3
2009 06 18 17 20 49	51.396	143.352	3.4	1.6
2009 10 15 09 52 47	51.276	143.022	8.3	2.2
2009 11 18 01 49 05	51.426	143.01	31.7	2.9
2009 11 18 06 57 05	51.364	143.161	13.9	1.8
2009 12 16 08 57 57	51.435	143.072	15.1	1.8
2009 12 16 11 29 45	51.433	143.088	18	1.7
2010 01 11 15 26 34	51.465	143.018	10	0.1
2010 08 23 13 36 11	51.33	143.05	3.5	0.8
2010 10 12 17 34 56	51.212	143.043	10	0.5
2011 03 15 10 53 07	51.58	143.212	14.3	2
2011 03 15 21 12 24	51.599	143.178	14.4	2.8
2011 04 26 22 05 08	51.29	143.021	15	1.3
2011 06 11 07 19 41	51.536	143.009	11.7	4
2011 07 21 17 22 57	51.468	143.124	15	2.2
2011 10 19 04 47 09	51.439	143.054	4	2.1
2011 11 10 12 37 04	51.531	143.189	14	1.8
2011 12 25 03 02 43	51.21	143.227	4.2	2.2
2012 01 27 19 02 37	51.42	143.194	17	1.3
2012 03 18 18 17 33	51.307	143.406	13	1.5
2012 05 28 10 26 26	51.396	143.092	9.3	1.8
2012 09 28 21 19 45	51.244	143.413	4.1	2.1
2012 10 05 03 55 07	51.325	143.272	8.3	1.5
2012 10 07 02 24 47	51.565	143.191	22.6	2.5
2012 10 07 03 03 13	51.593	143.173	24	2.2
2012 10 21 22 11 50	51.334	143.239	4.3	1.3
2012 12 28 06 53 11	51.351	143.091	5	2.5
2013 01 03 10 52 12	51.301	143.362	10	1.6
2013 01 16 13 41 36	51.221	143.017	4.7	1.7
2013 02 09 14 16 56	51.456	143.194	8	2.6
2013 02 11 04 18 22	51.471	143.207	5	2.9
2013 03 27 03 36 33	51.459	143.177	7.5	1.6
2013 04 05 14 32 25	51.506	143.073	15.1	1.6
2013 04 14 03 23 55	51.537	143.18	7.5	2.1
2013 05 18 21 23 00	51.133	143.368	15	2.4
2013 08 06 12 46 28	51.343	143.402	13.4	2.7
2013 08 06 15 44 24	51.309	143.404	4.2	1.1
2013 08 06 16 47 33	51.565	143.107	14.8	2.1
2013 08 09 17 46 33	51.566	143.102	14.6	1.5
2013 08 11 15 00 44	51.566	143.098	12.5	1.3
2013 09 20 03 10 35	51.404	143.459	9.4	3
2013 09 20 12 27 12	51.429	143.415	9.1	2.1
2013 09 20 21 16 10	51.424	143.397	5	1
2013 09 20 21 17 28	51.431	143.38	5	0.5
2013 10 14 02 31 04	51.424	143.441	9.8	3.3
2013 10 14 02 31 24	51.425	143.417	9.1	3.6
2013 10 14 02 32 02	51.424	143.448	9.3	3.2
2013 10 14 07 34 14	51.43	143.426	8.6	2.9
2013 10 14 08 36 24	51.42	143.463	8.9	4.3
2013 10 14 08 50 12	51.413	143.396	9.1	2
2013 10 14 16 09 42	51.431	143.439	8.4	2
2013 10 15 09 39 04	51.313	143.011	14.3	3.1
2013 10 15 16 23 17	51.313	143.032	13.6	1.8
2013 10 15 20 13 29	51.309	143.035	10.4	1.9
2013 10 16 02 12 48	51.411	143.359	9.8	2.2
2013 10 17 12 29 24	51.403	143.349	9.8	1.7
2013 12 19 14 30 38	51.101	143.66	10	2
2014 01 14 14 15 52	51.418	143.435	2	1.5
2014 01 20 13 52 32	51.235	143.269	4.1	1.2
2014 03 29 09 02 20	51.369	143.133	11.4	1.4
2014 06 12 07 50 40	51.564	143.109	18.7	1.4
2014 06 17 10 50 20	51.56	143.075	15.1	1.9
2014 06 17 14 30 32	51.428	143.401	9.3	1.6
2014 06 17 16 41 32	51.56	143.116	17.2	1.6
2014 06 25 13 19 57	51.4	143.122	15.1	1.3
2014 06 27 00 05 13	51.566	143.113	19.8	2.7
2014 06 27 15 38 46	51.567	143.106	18.7	3
2014 06 27 16 42 35	51.564	143.104	17.5	1.8
2014 06 27 20 41 27	51.567	143.103	18.2	1.7
2014 06 28 00 18 05	51.57	143.104	20.7	2
2014 06 29 16 04 43	51.566	143.096	17.3	2.1
2014 06 29 18 04 53	51.568	143.081	15	1.4
2014 06 29 22 46 45	51.569	143.097	20.9	2.1
2014 06 30 14 42 29	51.563	143.115	20.4	3.9
2014 06 30 14 48 12	51.563	143.108	19.3	2.6
2014 06 30 15 11 27	51.572	143.094	19.8	3.6
2014 06 30 16 19 49	51.572	143.095	18.2	2.3
2014 06 30 16 57 21	51.571	143.086	18.5	1.5
2014 06 30 17 40 39	51.577	143.096	17.8	3
2014 06 30 18 45 23	51.575	143.091	15	1.8
2014 06 30 19 28 22	51.541	143.005	3.9	1.4
2014 06 30 20 58 12	51.57	143.101	16.6	4.5
2014 06 30 21 06 12	51.565	143.086	19.7	2.9
2014 06 30 21 09 57	51.569	143.088	19.8	2.7
2014 06 30 21 18 52	51.504	143.095	21.9	1.8
2014 06 30 21 24 19	51.566	143.078	20	2.2
2014 06 30 21 25 08	51.577	143.08	18.8	2.7
2014 06 30 21 32 10	51.575	143.076	18.5	2.2
2014 06 30 21 37 43	51.575	143.083	19.9	2.6
2014 06 30 22 12 39	51.572	143.075	14.4	1.7
2014 06 30 22 23 45	51.571	143.074	20	2.3
2014 06 30 22 57 15	51.567	143.102	18.9	2.7
2014 06 30 22 57 15	51.571	143.1	20.5	2.7
2014 06 30 22 57 16	51.568	143.102	18.6	2.6
2014 06 30 23 25 13	51.562	143.107	19.1	2.7
2014 06 30 23 57 02	51.569	143.095	15.1	1.8
2014 07 01 00 00 24	51.565	143.102	18.9	1.7
2014 07 01 00 18 31	51.581	143.072	14.8	2
2014 07 01 00 48 38	51.561	143.1	18.8	3
2014 07 01 01 03 20	51.574	143.093	15.1	2.4
2014 07 01 01 14 11	51.579	143.081	13.5	3.2
2014 07 01 05 57 10	51.581	143.086	19.6	2.2
2014 07 01 05 57 10	51.579	143.081	18.7	2.2
2014 07 01 07 09 55	51.582	143.088	15	2.4
2014 07 01 07 42 07	51.574	143.076	10.1	1.5
2014 07 01 07 55 43	51.582	143.072	12.8	1.7
2014 07 01 09 20 37	51.579	143.074	18.3	2.1
2014 07 01 09 57 48	51.565	143.109	16.1	1.6
2014 07 01 10 21 27	51.563	143.083	18.8	1.9
2014 07 01 11 21 37	51.571	143.106	18.2	2.8
2014 07 01 12 30 49	51.565	143.084	18.5	2
2014 07 01 13 47 43	51.549	143.003	4.1	1.2
2014 07 01 14 32 18	51.577	143.092	17.8	2
2014 07 01 17 53 51	51.569	143.078	13.4	1.8
2014 07 01 18 03 06	51.568	143.079	17.3	1.8
2014 07 01 23 29 01	51.576	143.077	18	2.1
2014 07 02 00 16 15	51.566	143.105	17.3	2.4
2014 07 02 03 40 34	51.573	143.068	15	1.9
2014 07 02 10 53 33	51.569	143.074	19.1	2.4
2014 07 02 11 39 06	51.567	143.072	18.5	2.9
2014 07 02 16 13 19	51.561	143.096	18.4	1.4
2014 07 02 17 38 05	51.563	143.096	19.2	1.7
2014 07 02 19 02 47	51.562	143.1	19.3	2
2014 07 02 22 56 44	51.58	143.078	15	2.2
2014 07 03 00 26 54	51.586	143.093	19.2	4
2014 07 03 13 38 09	51.579	143.073	18.7	2.5
2014 07 03 14 28 43	51.56	143.108	18.6	1.8
2014 07 03 20 01 09	51.569	143.096	18.6	2.2
2014 07 03 20 29 56	51.564	143.109	18.2	2.7
2014 07 03 21 01 41	51.571	143.077	12.2	1.6
2014 07 03 22 42 17	51.583	143.104	13.5	1.9
2014 07 04 04 20 00	51.583	143.096	13.3	2.1
2014 07 04 04 27 34	51.585	143.1	14.4	1.8
2014 07 04 18 44 52	51.592	143.062	15	1.9
2014 07 04 21 10 58	51.583	143.1	14.1	1.9
2014 07 05 10 31 30	51.579	143.079	12.9	1.4
2014 07 05 14 27 36	51.589	143.077	15	2.7
2014 07 05 19 13 37	51.567	143.093	15.1	1.6
2014 07 06 08 45 08	51.6	143.085	17.1	1.8
2014 07 06 14 40 13	51.6	143.091	17.1	2.4
2014 07 07 01 20 09	51.561	143.075	15.1	1.9
2014 07 07 02 04 37	51.594	143.077	13.3	2.4
2014 07 07 13 48 04	51.59	143.079	15	1.8
2014 07 09 17 55 56	51.593	143.091	20.5	2.6
2014 07 09 20 13 12	51.594	143.065	15	2.3
2014 07 10 04 12 29	51.598	143.077	14.1	2.3
2014 07 10 12 26 27	51.599	143.089	18.2	1.9
2014 07 11 01 50 18	51.594	143.078	19.2	2.9
2014 07 11 03 38 51	51.592	143.068	11.3	1.9
2014 07 13 17 53 22	51.57	143.079	12.1	1.1
2014 07 14 08 44 55	51.589	143.085	18.2	2.1
2014 07 14 21 40 32	51.588	143.106	16.9	2.3
2014 07 14 21 42 05	51.588	143.104	15.9	2.1
2014 07 15 04 27 50	51.538	143.122	21.9	1.8
2014 07 15 10 04 35	51.587	143.1	13	1.8
2014 07 15 13 01 25	51.569	143.101	15	2.1
2014 07 15 19 58 33	51.561	143.103	16.8	1.4
2014 07 15 21 37 14	51.598	143.088	15	2.6
2014 07 15 22 40 57	51.596	143.094	16.3	1.9
2014 07 16 01 06 39	51.6	143.076	18.2	2
2014 07 16 09 50 17	51.595	143.087	14.1	2.1
2014 07 16 16 04 33	51.575	143.099	12.6	1.3
2014 07 17 08 19 37	51.595	143.087	17.9	2.1
2014 07 17 23 58 27	51.579	143.084	16.8	2.3
2014 07 18 01 14 00	51.577	143.069	12.7	1.5
2014 07 18 23 41 41	51.599	143.087	18.8	3.1
2014 07 19 15 08 56	51.597	143.092	15.1	2.8
2014 07 19 15 30 33	51.6	143.094	17.4	2.8
2014 07 19 15 30 33	51.595	143.096	15.4	2.7
2014 07 19 15 36 01	51.597	143.102	16.7	2.2
2014 07 19 15 45 05	51.596	143.089	15	1.8
2014 07 19 17 28 03	51.599	143.093	16.5	1.5
2014 07 19 17 29 54	51.597	143.097	17.3	2.3
2014 07 19 17 38 41	51.596	143.089	13.5	1.5
2014 07 22 05 31 19	51.564	143.093	12.3	2.2
2014 07 25 10 18 45	51.568	143.062	12.3	1.7
2014 08 02 22 45 36	51.594	143.082	18	1.8
2014 08 03 08 10 10	51.58	143.072	15.7	1.8
2014 08 03 10 38 38	51.585	143.077	17.7	2
2014 08 03 10 57 47	51.582	143.061	22.1	1.4
2014 08 03 16 04 11	51.558	143.11	17.3	1.4
2014 08 06 18 40 31	51.567	143.1	13.6	1.3
2014 08 11 18 03 35	51.584	143.081	10.5	1.1
2014 08 17 17 31 08	51.569	143.09	17.2	1.6
2014 08 21 13 17 52	51.299	143.404	15	1.8
2014 09 04 12 56 11	51.546	143.189	15	1.7
2014 09 14 15 22 38	51.551	143.17	15	2
2014 09 14 15 25 00	51.543	143.199	15	1.6
2014 09 16 15 13 21	51.572	143.192	15.1	1.5
2014 09 19 08 16 42	51.351	143.112	10.6	1.7
2014 10 08 00 54 35	51.587	143.102	13.2	1.8
2014 10 08 02 39 18	51.561	143.091	12.1	2
2014 10 30 16 10 36	51.408	143.5	5.4	2
2014 10 31 12 13 28	51.582	143.11	12.8	1.6
2014 10 31 23 26 03	51.581	143.113	12.9	2
2014 11 01 16 59 11	51.573	143.107	13.3	2.4
2014 11 02 12 50 41	51.551	143.04	4.1	1.3
2014 11 06 20 05 57	51.57	143.124	13.1	1.7
2014 11 12 18 21 16	51.257	143.49	4.1	2
2014 12 10 20 54 22	51.542	143.089	2.4	1.3
2014 12 21 10 59 41	51.481	143.601	12.3	2.3
2014 12 22 17 13 40	51.494	143.622	7.6	2.3
2014 12 29 19 06 07	51.488	143.607	12.2	1.7
2015 01 01 14 23 57	51.372	143.137	4	1.1
2015 01 09 09 43 24	51.567	143.091	12.6	1.6
2015 01 14 08 17 40	51.519	143.599	7.5	1.8
2015 01 18 05 55 30	51.578	143.095	11	1.3
2015 01 18 06 05 47	51.578	143.094	11	1.7
2015 01 18 08 54 22	51.58	143.097	11	1.4
2015 01 18 08 54 22	51.577	143.094	11.1	1.4
2015 01 18 15 52 37	51.578	143.096	14.4	2.2
2015 01 19 12 43 24	51.417	143.118	12.3	1.2
2015 01 21 13 59 41	51.181	143.103	4.2	1.2
2015 01 29 13 57 52	51.575	143.086	12.8	1.3
2015 01 31 07 04 21	51.6	143.095	15.8	1.2
2015 02 21 04 44 01	51.58	143.083	12.4	1.4
2015 03 07 03 07 24	51.427	143.117	13.3	1.3
2015 03 07 20 07 57	51.41	143.334	14.2	0.8
2015 03 14 23 49 11	51.564	143.131	4.4	1.8
2015 04 04 10 44 30	51.56	143.112	13.2	2.1
2015 04 04 20 00 48	51.562	143.113	12	1.1
2015 04 04 20 14 16	51.561	143.111	11.1	1.1
2015 04 05 18 46 20	51.558	143.11	11.9	1.8
2015 04 08 15 11 38	51.563	143.085	11.4	1
2015 05 14 16 44 23	51.445	143.066	11.3	1.2
2015 06 29 17 34 29	51.343	143.487	10.4	2
2015 06 29 17 55 32	51.356	143.506	9.7	1.7
2015 07 15 13 57 07	51.166	143.004	5.9	1.5
2015 08 07 13 57 02	51.591	143.169	4.1	1.2
2015 08 08 18 45 18	51.552	143.181	15.1	1.8
2015 08 08 18 50 02	51.545	143.195	15.1	2.6
2015 08 08 21 53 21	51.551	143.176	15	1.9
2015 08 14 12 02 30	51.556	143.204	15.1	1.8
2015 08 14 19 39 41	51.597	143.151	7.3	1.3
2015 08 21 02 34 35	51.578	143.204	14.4	2.4
2015 08 30 04 03 07	51.54	143.216	15	2.1
2015 08 30 04 38 21	51.558	143.175	4.2	1.6
2015 08 30 16 54 51	51.212	143.091	10	1.6
2015 09 09 11 56 53	51.539	143.269	12	1.3
2015 09 23 09 06 58	51.429	143.108	15	1.9
2015 09 26 00 54 31	51.395	143.13	4.3	1.3
2015 10 08 22 46 56	51.588	143.137	12.2	1.9
2015 10 09 01 32 04	51.583	143.148	10.6	2.1
2015 11 24 18 22 25	51.506	143.473	4.1	1.5
2015 11 27 16 30 30	51.462	143.081	16.2	1.7
2015 12 18 07 49 06	51.392	143.345	27.6	1.6
2015 12 31 04 17 08	51.524	143.171	10.5	2.8
2015 12 31 19 42 27	51.265	143.022	15.1	1
2016 01 02 21 43 14	51.53	143.183	4.3	0.7
2016 01 13 09 25 05	51.417	143.155	15.9	1.5
2016 01 13 10 05 29	51.299	143.078	4.1	1.1
2016 01 13 16 22 07	51.171	143.002	5	1.2
2016 01 15 20 03 22	51.337	143.498	12	1.3
2016 02 07 23 59 52	51.353	143.42	11.3	2
2016 02 29 21 31 25	51.488	143.081	4.2	1.1
2016 03 01 09 55 41	51.544	143.076	12.3	1.1
2016 03 01 12 53 18	51.595	143.111	17	1.8
2016 03 02 08 43 51	51.316	143.234	5	1
2016 03 05 09 32 45	51.489	143.017	3.3	2.7
2016 03 05 11 30 51	51.498	143.006	4	0.8
2016 03 05 11 46 22	51.5	143.008	8	1.2
2016 03 18 09 59 06	51.301	143.272	5	0.7
2016 03 19 00 51 05	51.472	143.088	12.5	2.1
2016 03 24 19 06 48	51.55	143.091	9	2.3
2016 03 24 23 58 25	51.537	143.074	0.6	1.3
2016 04 11 19 04 41	51.369	143.042	10	1.6
2016 04 18 18 54 45	51.565	143.086	15.2	1.1
2016 04 20 21 50 09	51.245	143.338	5	1.5
2016 04 22 05 49 26	51.224	143.181	2	1.9
2016 04 25 18 45 17	51.241	143.081	10	1.7
2016 04 29 15 26 03	51.492	143.608	10	1.1
2016 05 23 12 57 16	51.57	143.136	2.8	1.7
2016 07 20 14 57 15	51.314	143.187	10	1.8
2016 07 29 13 36 29	51.309	143.355	5	1.3
2016 07 30 17 35 14	51.4	143.16	15.8	1.3
2016 08 01 17 28 39	51.446	143.437	5	0.4
2016 08 01 19 25 56	51.441	143.415	6.5	3.1
2016 08 01 19 33 40	51.439	143.424	6.8	1.7
2016 08 03 09 44 56	51.459	143.428	5	0.8
2016 08 04 02 06 13	51.292	143.272	5	1.1
2016 08 05 08 08 25	51.411	143.353	4.1	1.4
2016 08 08 12 49 39	51.426	143.349	10	1.3
2016 08 11 17 34 20	51.527	143.053	6.8	1
2016 08 31 17 45 17	51.459	143.032	10	1.1
2016 10 07 20 25 53	51.344	143.055	11.3	1.7
2016 10 13 22 20 14	51.566	143.134	10.2	1.6
2016 10 26 13 00 49	51.377	143.071	13.5	1.3
2016 11 26 10 31 43	51.362	143.011	16	1
2016 11 27 08 59 17	51.523	143.201	4.3	1.3
2016 11 29 00 08 04	51.6	143.097	13.7	1.8
2016 12 11 07 38 02	51.123	143.201	10	1.5
2016 12 11 17 50 20	51.416	143.596	9.6	1.6
2016 12 20 13 27 53	51.294	143.54	13.5	2.4
2016 12 24 00 30 03	51.557	143.121	11.9	1.7
2016 12 26 22 40 38	51.369	143.043	11.6	1.3
2017 01 09 12 09 36	51.372	143.279	14.5	1
2017 02 21 09 24 04	51.573	143.123	2.4	1.1
2017 02 25 16 27 24	51.51	143.457	21.3	1.2
2017 02 27 14 27 39	51.354	143.031	15	1.1
2017 04 30 08 56 14	51.467	143.016	19	2.8
2017 04 30 08 56 30	51.519	143.039	19	2.6
2017 05 09 23 15 32	51.592	143.315	4.2	2.5
2017 06 10 14 38 37	51.401	143.521	5.5	1.4
2017 07 08 15 39 27	51.245	143.526	12	2.7
2017 07 12 10 57 52	51.523	143.055	1.2	1.1
2017 07 29 22 53 47	51.409	143.371	3.9	1.1
2017 08 28 13 57 55	51.311	143.534	10.5	1.7
2017 08 31 18 42 21	51.409	143.695	10.4	2.3
2017 09 30 00 49 59	51.35	143.408	8.4	1.4
2018 01 21 23 48 06	51.524	143.086	11	1.2
2018 01 22 01 17 18	51.513	143.088	1.5	1.3
2018 02 09 19 37 01	51.338	143.046	4.3	1.1
2018 04 06 20 42 10	51.515	143.535	0.5	2
2018 04 09 16 06 16	51.376	143.329	5	0.9
2018 04 19 11 36 27	51.573	143.133	2.7	0.7
2018 04 24 01 16 45	51.557	143.049	4.3	1.3
2018 04 26 02 29 55	51.559	143.114	2.6	1.4
2018 04 29 04 27 18	51.447	143.548	13	1.9
2018 05 24 17 41 44	51.581	143.079	10.6	2.1
2018 05 28 11 50 26	51.57	143.055	4.2	1.3
2018 06 21 15 33 07	51.468	143.475	4.1	1.4
2018 07 05 22 49 44	51.524	143.187	5	1.8
2018 07 11 12 07 16	51.368	143.489	2.1	1.2
2018 07 15 19 40 42	51.475	143.044	4.2	0.9
2018 07 17 19 40 53	51.36	143.111	8	1.4
2018 07 23 12 54 00	51.254	143.649	10.9	1.8
2018 07 31 10 05 59	51.345	143.267	20	1.6
2018 08 19 11 48 06	51.379	143.537	11.3	1.7
2018 08 28 14 53 17	51.168	143.012	5	1.3
2018 09 11 08 45 33	51.347	143.431	0	1.5
2018 09 18 08 20 40	51.336	143.046	8.5	1.9
2018 09 19 07 55 26	51.336	143.053	4.2	1.4
2018 09 20 07 36 24	51.199	143.201	6.6	1.5
2018 09 25 12 31 14	51.33	143.045	10.4	1.1
2018 10 05 22 01 08	51.375	143.531	11.5	1.7
2018 10 20 14 02 30	51.463	143.055	4.2	1.2
2018 10 28 22 48 13	51.513	143.011	4.2	1.6
2018 11 04 14 47 22	51.381	143.543	9.9	1.8
2018 11 07 19 47 08	51.544	143.178	4.3	1.3
2018 12 03 12 44 10	51.565	143.108	16.2	1.2
2018 12 12 04 59 39	51.223	143.197	2	1.6
2018 12 20 12 17 47	51.374	143.521	10.1	1.3
2018 12 25 15 02 52	51.486	143.475	11.7	3.3
2018 12 25 15 41 01	51.487	143.449	11	3.2
2018 12 30 15 28 31	51.423	143.109	11.9	2
2019 02 19 17 54 51	51.38	143.525	11	0.8
2019 02 27 22 54 45	51.311	143.406	9.6	2.3
2019 03 10 14 46 54	51.456	143.433	14.1	1.4
2019 03 12 15 31 58	51.171	143.937	14.8	1.4
2019 03 19 13 31 31	51.317	143.021	4.2	0.7
2019 04 01 17 31 39	51.563	143.091	2.8	1.1
2019 04 09 19 50 10	51.221	143.056	1.2	1.1
2019 05 01 01 25 30	51.369	143.145	11.8	1.5
2019 05 02 04 13 33	51.589	143.101	15	1.3
2019 05 02 04 43 57	51.232	143.483	0.6	2.1
2019 05 08 13 54 40	51.414	143.147	15.7	1.4
2019 05 13 15 58 32	51.518	143.468	11.8	2.3
2019 05 24 15 48 11	51.268	143.553	7.7	1.5
2019 06 06 23 58 59	51.498	143.527	3.8	1.7
2019 06 17 01 32 22	51.398	143.126	12.1	1
2019 07 03 15 45 47	51.454	143.04	5	0.5
2019 08 05 14 08 03	51.401	143.172	4.1	1.5
2019 08 31 16 32 52	51.468	143.442	0	1
2019 09 03 11 41 07	51.57	143.131	17.5	2.1
2019 09 12 14 20 22	51.458	143.544	5	1.1
2019 10 15 10 02 00	51.568	143.12	12.2	1.9
2019 11 24 16 06 43	51.18	143.155	2.8	1.5
2019 12 05 03 57 08	51.553	143.007	27.6	1.4
2019 12 30 18 22 33	51.554	143.104	11.7	2.6
2019 12 30 18 24 27	51.557	143.102	15.8	1.3
2020 01 01 12 54 47	51.561	143.103	13.1	0.8
2020 01 05 09 39 52	51.563	143.097	14.8	0.9
2020 01 05 13 05 54	51.559	143.091	15	1.1
2020 01 05 13 09 14	51.564	143.1	14.6	0.8
2020 01 05 14 40 37	51.557	143.099	11.4	2.2
2020 01 05 16 14 50	51.562	143.097	12.4	2.3
2020 01 05 16 16 20	51.562	143.105	15	1.1
2020 01 05 16 41 29	51.56	143.101	15	2.6
2020 01 05 16 52 01	51.562	143.1	18.4	1.1
2020 01 06 13 45 28	51.566	143.088	13.8	1.1
2020 01 06 16 38 21	51.565	143.079	16	1
2020 01 06 17 36 38	51.562	143.076	16	0.8
2020 01 09 17 43 21	51.562	143.095	10	1.2
2020 01 11 04 25 35	51.505	143.06	16	1.3
2020 01 12 06 09 02	51.192	143.202	4.1	1.4
2020 01 25 14 43 25	51.369	143.104	11.9	0.9
2020 01 30 14 11 26	51.363	143.107	9.9	1.2
2020 02 05 22 31 07	51.572	143.086	2.3	1.7
2020 02 25 18 05 34	51.323	143.553	1.1	1.4
2020 02 25 18 39 16	51.325	143.547	1.2	1.4
2020 03 02 18 13 25	51.317	143.56	5	1.1
2020 03 22 14 45 13	51.595	143.141	15	1.6
2020 03 28 20 59 14	51.393	143.181	12.9	1.9
2020 03 30 07 30 20	51.398	143.183	7.5	2.6
2020 04 27 20 51 13	51.581	143.081	11.4	2.1
2020 05 13 17 58 42	51.56	143.103	17.1	2.7
2020 05 14 19 02 06	51.561	143.096	15.7	1.8
2020 05 28 12 50 16	51.561	143.093	14.2	1.2

magnitude_det_rate_estimation_mle.py

"""
File: magnitude_det_rate_estimation_mle.py
Author: Andrey Stepnov
Email: a.stepnov@geophystech.ru
Github: https://github.com/jamm1985

Description: Estimation the detection bound of local magnitude using a widely used statistical model developed by Ogata and Katsura. **Constant** rate is assumed.

Ogata, Y., Katsura, K. (1993). https://doi.org/10.1111/j.1365-246x.1993.tb04663.x

Input: Earthquake catalog produced by REPORT command from SEISAN distribution 

https://seis.geus.net/software/seisan/node124.html

Output: Estimated mu (50% detection rate), b-value and sigma. Fig_1.png with the magnitude-frequency relations: ‘+’ are observed magnitudes, line is estimated frequency-densities curves for catalog. X-axes is for local magnitude and Y-axes given in logarithmic scale.

Dependencies: scipy >= 1.7.3 
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
from scipy import special
from scipy.stats import norm
from scipy.optimize import minimize
import argparse

parser = argparse.ArgumentParser(description='Process input')
parser.add_argument('report', type=str, help='Input report.out')
parser.add_argument(
    '--min_mag',
    type=float,
    default=1.0,
    help='min_mag for b-value estimation (default: 1.0)'
)

args = parser.parse_args()

input_report = args.report
min_mag = args.min_mag

def density(M,theta):
    """evaluation function from observation and given parameters
    :M: M_1,...,M_n
    :theta: parameter vector (mu, betta, and sigma)
    :returns: probability for every observation
    """
    beta = theta[0]
    mu = theta[1]
    sigma = theta[2]
    constant = (np.exp(-beta*M)*norm.cdf(M,mu,sigma)).sum()
    return np.exp(-beta*M)*norm.cdf(M,mu,sigma)/constant


def log_l(theta):
    """log-likelihood function from observation and given parameters
    :M_observation: M_1,...,M_n (global context!)
    :theta: parameter vector (mu, beta, and sigma)
    :returns: log-likelihood value for all observations
    """
    beta = theta[0]
    mu = theta[1]
    sigma = theta[2]
    N = M_observation.size
    return (N * np.log(beta) - beta * (M_observation + np.log(norm.cdf(M_observation,mu,sigma))).sum() + np.exp(N*beta*mu - N/2 * beta**2 * mu**2))


def der_log_1(theta):
    """first derivatives of log_l
    :M_observation: M_1,...,M_n (global context!)
    :theta: parameter vector (mu, beta, and sigma)
    :returns: gradient vector with values of first derivatives
    """
    beta = theta[0]
    mu = theta[1]
    sigma = theta[2]
    N = M_observation.size
    der = np.zeros_like(theta)
    exp_term = np.exp(N*beta*mu - N/2 * beta**2 * mu**2)
    exp_term_beta_prime = N*mu - N*beta*mu**2
    exp_term_mu_prime = N*beta - N*beta**2*mu
    der[0] =  N/beta - (M_observation + np.log(norm.cdf(M_observation,mu,sigma))).sum() + exp_term*exp_term_beta_prime
    der[1] =  N*beta*(1/sigma) + exp_term * exp_term_mu_prime
    der[2] =  beta* (((M_observation - mu)/sigma**2).sum())
    return der

def hess_log_1(theta):
    """second derivatives of log_l
    :M_observation: M_1,...,M_n (global context!)
    :theta: parameter vector (mu, beta, and sigma)
    :returns: hessian matrix with values of second derivatives
    """
    beta = theta[0]
    mu = theta[1]
    sigma = theta[2]
    N = M_observation.size
    hess = np.zeros((3,3))
    exp_term = np.exp(N*beta*mu - N/2 * beta**2 * mu**2)
    exp_term_beta_prime = N*mu - N*beta*mu**2
    exp_term_mu_prime = N*beta - N*beta**2*mu

    hess[0,0] = -N/beta**2 + exp_term_beta_prime**2*exp_term - N*mu**2*exp_term
    hess[0,1] = N*(1/sigma) + exp_term_beta_prime*exp_term*exp_term_mu_prime + exp_term*(N-2*N*beta*mu)
    hess[0,2] = ((M_observation - mu)/sigma**2).sum()

    hess[1,0] = N/sigma**2 + exp_term_mu_prime*exp_term*exp_term_beta_prime + exp_term*(N-2*N*beta*mu)
    hess[1,1] = exp_term_mu_prime**2*exp_term - N*beta**2*exp_term
    hess[1,2] = -beta*N*(1/sigma**2)

    hess[2,0] = ((M_observation - mu)/sigma).sum()
    hess[2,1] = -N*beta*(1/sigma)
    hess[2,2] = -beta*((M_observation - mu)/sigma**2).sum()

    return hess

# estimate b-value 
# this is modified version of 
# source: https://github.com/qingkaikong/Learning_Obspy/blob/master/04_catalogAnalysis_20140421/calc_rec.py
def calc_recurrence(magnitudes, min_mag = None, max_mag = None, interval = 0.05):

    """
    This function reads an earthquake catalogue file and calculates the
    Gutenberg-Richter recurrence parameters using maximum likelihood (Aki 1965) 
    approach.

    Funtion arguments:

     magnitudes: list of magnitudes values
        min_mag: minimum magnitude for which data will be used - i.e. catalogue
                completeness
        max_mag: maximum magnitude used in bounded G-R curve. If not specified,
                defined as the maximum magnitude in the catlogue + 0.1 magnitude
                units.
        interval: Width of magnitude bins for generating cumulative histogram
                of earthquake recurrence. Default avlue 0.1 magnitude units.  
    
    """
    
    # If minimum magnitude is not specified, read all magnitudes
    if min_mag is not None:
        pass
    else:
        min_mag = -1.0
        
    # If minimum magnitude is not specified default value to minimum in catalogue
    if min_mag == -1.0:
        min_mag = min(magnitudes)
    # If maximum magnitude is not specified default value to maximum in catalogue
    if max_mag is not None:
        pass
    else:
        max_mag = max(magnitudes) + 0.1

    num_eq = len(magnitudes)
    print('Minimum magnitude: {}'.format(min_mag))
    print('Total number of earthquakes: {}'.format(num_eq))
    #num_years = max(years)-min(years)
    #annual_num_eq = num_eq/num_years
    #print 'Annual number of earthquakes greater than Mw', min_mag,':', \
    #annual_num_eq
    print('Maximum catalog magnitude: {}'.format(max(magnitudes)))
    print('Mmax = {}'.format(max_mag))
    max_mag_bin = max(magnitudes) + 0.15
    
    # Magnitude bins
    bins = np.arange(min_mag, max_mag_bin, interval)
    # Magnitude bins for plotting - we will re-arrange bins later
    plot_bins = np.arange(min_mag, max_mag, interval)

    ###########################################################################
    # Generate distribution
    ###########################################################################
    # Generate histogram
    hist = np.histogram(magnitudes,bins=bins)

    # Reverse array order
    hist = hist[0][::-1]
    bins = bins[::-1]

    # Calculate cumulative sum
    cum_hist = hist.cumsum()
    # Ensure bins have the same length has the cumulative histogram.
    # Remove the upper bound for the highest interval.
    bins = bins[1:]

    # Maximum Likelihood Estimator fitting
    # b value
    b_mle = np.log10(np.exp(1)) / (np.mean(magnitudes) - min_mag)
    beta_mle = np.log(10) * b_mle
    print('Maximum Likelihood: b value {}'.format( b_mle))

    return b_mle, cum_hist[::-1], bins[::-1]


# get M_1,...,M_n from catalog
catalog = pd.read_fwf(input_report)
M_observation = catalog['Ml'].dropna().to_numpy()

# obtain initial theta0
mu_0 = M_observation.mean()
b_value_0, _, _ = calc_recurrence(M_observation.tolist(), min_mag = min_mag, max_mag = None, interval = 0.05)
beta_0 = np.log(10)*b_value_0
std_0 = M_observation.std()
theta0=[beta_0,mu_0,std_0]
bounds=((0.3*np.log(10),0.8*np.log(10)),(3.0,5.0),(0.01,0.8))
print('theta0 beta, mu, sigma is {}'.format(theta0))
print('bounds is {}'.format(bounds))

# estimate numerically using trust region method
res_exact = minimize(log_l, theta0, method='trust-exact', jac=der_log_1, hess=hess_log_1, options={'gtol': 1e-8, 'maxiter': 10000, 'disp': True}, bounds=bounds)
print('estimated theta is {}'.format(res_exact.x))

# estimate using constrains
res_constrains = minimize(log_l, theta0, method='trust-constr', jac=der_log_1, hess=hess_log_1, options={'xtol': 1e-8, 'verbose': 1, 'maxiter': 10000}, bounds=bounds)
print('Current function value:{}'.format(log_l(res_constrains.x)))
print('estimated theta is {}'.format(res_constrains.x))

# choose final theta
if log_l(res_exact.x) < log_l(res_constrains.x):
    final_theta = res_exact.x
    print('exact method is better, final theta is {}'.format(final_theta))
else: 
    final_theta = res_constrains.x
    print('bounded method is better, final theta is {}'.format(final_theta))

# plot results
plt.cla()
plt.rcParams.update(plt.rcParamsDefault)
plt.rcParams['figure.figsize'] = [10, 6]
plt.rcParams['ytick.major.size'] = 10
#plt.rcParams['ytick.major.width'] = 4
plt.rcParams['ytick.minor.size'] = 5
#plt.rcParams['xtick.minor.width'] = 2
#plt.title("Prior 2006-09",fontsize = 15)
plt.ylabel('Probability', fontsize = 15)
plt.xlabel('Magnitude', fontsize = 15)
plt.yscale('log')
plt.yticks(fontsize=14)
plt.xticks(fontsize=15)
plt.ylim([0.001, 1])
plt.xlim([M_observation.min()-0.5, M_observation.max()+0.5])
plt.text(0.5, 0.5,r'$b$={}'.format(np.round(final_theta[0]/np.log(10),decimals=1))+'\n'+r'$\mu$={}'.format(np.round(final_theta[1],decimals=1))+'\n'+r'$\sigma$={}'.format(np.round(final_theta[2],decimals=1)), ha='left', va='center', fontsize = 15)
# ECDF
x = np.sort(M_observation)
y = np.arange(len(x))/float(len(x))
# EPMF
y[1:] -= y[:-1].copy()
u, inv = np.unique(x, return_inverse=True)
sums = np.zeros(len(u), dtype=y.dtype) 
np.add.at(sums, inv, y)
plt.scatter(u, sums, marker='+')
# estimated density
#x3=np.linspace(0,5,len(x))
plt.plot(u,density(u,final_theta))
plt.savefig('fig_1.png',dpi=600)