Statistical models for seismic magnitude

Physics of The Earth and Planetary Interiors, 21(1980) 237—260 © Elsevier Scientific Publishing Company, Amsterdam — Printed in The Netherlands

237

STATISTICAL MODELS FOR SEISMIC MAGNITUDE ANDERS CHRISTOFFERSSON Department of Statistics, University of Uppsala, Uppsala (Sweden)

(Received October 31, 1978; revised and accepted April 27, 1979)

Christoffersson, A., 1980. Statistical models for seismic magnitude. Phys. Earth Planet. Inter., 21: 237—260. In this paper some statistical models in connection with seismic magnitude are presented. Two main situations are treated. The first deals with the estimation ofmagnitude for an event, using a fixed network of stations and taking into account the detection and bias properties of the individual stations. The second treats the problem of estimating seismicity, and detection and bias properties of individual stations. The models are applied to analyze the magnitude bias effects for an earthquake aftershock sequence from Japan, as recorded by a hypothetical network of 15 stations. It is found that network magnitudes computed by the conventional averaging technique are considerably biased, and that a maximum likelihood approach using instantaneous noise-level estimates for non-detecting stations gives the most consistent magnitude estimates. Finally, the models are applied to evaluate the detection characteristics and associated seismicity as recorded by three VELA arrays: UBO (Uinta Basin), TFO (Tonto Forest) and WMO (Wichita Mountains).

1. Introduction One of the main problems in the estimation of magnitudes of seismic events from network data is that these events are not always detected by all stations. The standard procedure for estimating magnitude, by averaging the observed magnitudes of recording stations, gives estimates that are biased upwards. Methods to cope with this problem have been given by Herrin and Tucker (1972) who computed the expected error introduced by the averaging procedure, and Ringdal (1976) who developed a maximum likelihood procedure. The maximum likelihood method given in this paper differs from Ringdal’s approach, in that it takes into account the probability that the event is detected by the network. This will have an effect on the estimated magnitude for small events (events that are detected by a small number of stations). The maximum likelihood estimator for magnitude requires knowledge of the station detection parameters (threshold and slope of detection curve) and region—station bias.

In Section 3, methods to estimate these parameters are developed. This is done assuming that the magnitude distribution in the source region can be approxi. mated with the usual linear relationship between magnitude and logarithmic frequency. The first method uses single-station data for simultaneous estimation of the slope of the seismicity curve and the detection parameters of the station. This maximum likelihood estimator was, although not explicitly stated, given by Kelly and Lacoss (1969). One interesting result is that the distribution of log(A/T), rather than magnitude at a station, is independent of the focal locations of the events and scattering. The second method is developed for estimation of region—station bias and scattering standard deviation. Because of the complexity of the joint distribution in the general case, detailed calculations are carried out for only two stations. By combining analyses made for different combinations of station pairs, it is possible to estimate relative bias and scattering standard deviations for all stations in a network.

238

2. Theory

station. Then, we define the variable a1 such that

2.1. Basic assumptions It is assumed that the amplitude of the signal generated by a seismic event in a specific region and arriving at a station can be modelled as

a1 =Yj if the event is seen at station i, and a1 EN, if the event is not seen at station i. Unless otherwise stated, the distance—depth corrections (Q1) are regarded as known constants. The conditional distribution of a/rn is given by

A/T

h1~a,/m)da,

(1)

emeQeBee

where AlT is amplitude over period, rn is the “true” magnitude of the event, Q is the distance—depth correction, B is the region—station bias, i.e. the average scattering effect due to inhomogeneities in the Earth, and is the difference between total scatteringand average scattering, Taking logarithms of both sides in (1), we get log(A/T) = rn + Q i-B + e or, letting y we get y

=

=

(2)

=

[

(fi~ai/rn~~~ai G~)/’y~]da,for a1 ~N, (B, + Q, + m G,) —

—

_______

+~ If the scattering effects, ~, the joint distribution of (a1, a2, ~,

(3)

If we regard e as the sum of effects from many travel path inhomogeneities or scattering sources it can be shown (under general assumption) that the distribution of eis Gaussian. The conditional distribution of y for given “true” magnitude m is then 2)[y (B + Q + m)]2} lily/rn) = (l/~/~a) exp{—(l/2a (4) Here Q is assumed to be a known constant, and a2 is

M

h(a1,a2,

aMIm)dalda2

...,

F(detect/m) = 1

n ~ (—(B, rn ‘~ ~Ja~+7~ M

—

L

+

i=1

+

—

Gi))

M

=

1

—

Ii

F(a, E N,/rn)

(8)

j=1

The distribution of ~, au/rn) given that the event is detected by the network is then H(a 1, a2, aMIm)dal da~ ...

h(a1, 1

...

aM/m)dal ... da~ Al P(a, E Nj/rn)

...,

= —

y

=

LI h,(a,/m)da, (7)

If an event is declared as detected by the network whenever it is detected by at least one station, we obtam the detection probability for the network as

...,

—

da~=

...

/=1

—

the variance of the scattering, i.e., of e. We wifi further assume that the station has a detection curve giving the probability that the event is seen given y. We will here assume that this curve is of the form

for a~EN,

CJL,f are independent, ..., aM/rn) is

...,

log(A/T) denote the station magnitude,

rn + B + 12 + C

$ [Cv G)/’i]

]

(6)

n

(9)

1=1

(1 ~

~

2/272]dz (5) X exp [—(zbe—interpreted G) where G can as the average threshold and 7 as the standard deviation of the threshold. 2.2. Distribution of observed log(A/T) at a network

This is the conditional distribution of log(A/T) recorded at magnitudes the network.is Ifthe distribution of true earthquake known, we can derive the unconditional distribution of log(A/T) recorded at the network. Let v(rn)drn denote the distribution of earthquake magnitude in a region. The unconditional distribution of a 1, aM is ...,

of M stations

g(a1, a2, Let N, denote the subset of station log(A/T) for given true magnitudes that are not seen at the ith

=

...~

aM)dal

f h(a1,

...,

...

da~

aMIm)dal

...

daMv(m)drn

(10)

239

The probability that all a, EN,, i.e. the probability of not recording is F(a1 EN1,

a~jE NM)

...,

M

_

iT! ~ i=1

=

—(rn + B + Q,

—

~/ai + ~

G) v(m)din

(11)

giving the unconditional distribution of log amplitude recorded at the network G(ai, —

with 17(m,/m) f(rn1 + Qj/m). From (16) we see that the distribution of mg/rn depends on the isdistance—depth correction. However, this correction introduced in order to obtain mea-

aM)dal da~ g(a1, aM)dal da~q 1 —F(a1 EN1, ...,aMENM) ...~

...

...,

...

(12)

It is often assumed that there is a linear relation between magnitude and logarithmic frequency of earthquake occurrence. This corresponds to the assumption that the distribution of earthquake magnitude is u(rn)dm = ~3exp[~3(rn

—

rn0)] drn form> m0

u(m)drn = 0

form
(13)

sures of event magnitudes that are independent of the focal locations. In order to avoid this dependence, we can regard the distance—depth correction as a random variable, normally distributed with expectation ~ and standard deviation a~.We can then obtain the probability of detecting an event for a given station magnitude. This curve takes the same form as (5). We get the detection curve ~[(rn,G;)/7fl withG G, — Q,,’y’ =(y~~~Q) 2 ‘1/2

(17) and the distribution of m,/rn as h,’~(rn,/m) =

ff,~’(mi/rn)c~[(m, G)/7*Jdrnj form, ~N, i-rn G)/Va? + ~ for m1EN1 —

—

Inserting this in (12) we get

(18)

G(aI,...,aM)dal...daM

f h(a~, a0/m)da ... ,

...

Using h* instead of h, in (9) and (14) then gives the conditional and unconditional distribution of observed station magnitudes.

da~exp(—~3rn)dm (14)

mO

M

f (1

2.3. Comments —

11a~EN,/m) exp(_~n))dm

m0

The integrals in (14) exist for all m0 and are convergent when m0 -~ In many cases it may be desirable to consider the distribution of station magnitudes instead of station log(A/fl. Letting rn denote the observed magnitude at the ith station, we have —°°.

m, = a,

—

(15)

Q,

and the distribution of rn,/rn is then 4’(rn,/rn)dmi= hg(m, + Q’/rn)dm, h~ ~

=

~( ¶J~+

—(B,+Q,+rn—Gj)\ 7~

)

We have here defmed detection of an event as seen by at least one station. The use of other possible definitions, for example, events seen by at least two stations,will only change the denominators in (9) and (14), which in turn are derived from eq. 8. In the following sections most results are derived using log(A/T) as a basis. The corresponding results in terms of magnitude are obtained by substituting magnitude for log(A/7) and setting all distance—depthcorrections equal to zero. The station parameters G and 7 are then interpreted as detection parameters in terms of magnitude, andsimilar can betotransformed in a way that in eq. back 17. to a log(A/T) basis

form,~N,

form,EN 1 (16)

3. Estimation The two distributions (9) and (14) are in general suited for two different estimation problems. The con-

240

ditional case is useful mainly for estimation of magnitudes of individual events, while estimation of structural parameters such as seismicity, station bias, etc., is most conveniently done in the framework of the unconditional approach. In certain cases, the conditional

Frequency 20

1 5

approach used to estimate eters, and can for be those cases it has thestructural advantageparamof not

1~

requiring with Here theboth conditionaL any approaches knowledge of arelikelihood prior considered, distributions. beginning 3.1. Conditional maximum Consider first the case of one station and one event. In this case, the distribution of observed log(A/7) = a for a given true magnitude m ~ H(a/m)da

=

~ /~(B+ Q + rn

j~a/m)

~

/

~

~2

=

(19)

a) exp {— [a — (B + Q + rn)] 2/202 }

(1 /s./

f ((1/\/~) f

exp(-t2/2)dt)

________ —

—

zero shows It is shown probability that the in Appendix effect of notofdetecting non-detectionis 1 that an event.substantial. Figure 1 E(a/m) = B + Q + rn + (a2R,/~T

2)Z(x)

readily shown that the maximum likelihood estimate corresponds to solving the following equation for rn

log (A/I) 10

05

~

/ -05

/

x

f

exp(_x2/2)/~~=

(22)

In the case of a network of several stations, one might estimate the magnitude for each station that detected the event, according to the above approach, and then average these values to obtain an estimate of the true magnitude. However, if the stations that did

(21)

7

withxas before and

~

Fig. 1. Frequency distribution H(a/m) for m = 3.8, G = 0.0, B + Q = —4.0, a = 0.3, 7 = 0.2 and the corresponding Gaussian distribution with expectation —0.2 and standard deviation 0.3.

(20)

withy (a G)/7andx (B + Q i-rn G)/~/~+ 72~ The shape of the distribution is shown in Fig. 1, together with the corresponding Gaussian distribution, which corresponds to the case where the station has a

Z(x) =

05

a = B + Q + m + (a2/~,/~~)Z(x)

X (1 /~.J~)exp(—t2/2)dt /

6

-

+ 72

or explicitly

H(a/m)da

Gou AaIml -1 -05

0

0

exp(—t2/2)dt

/

-10

/

/

/

/

I I I I

/ -

Figure 2 shows E(a/rn) for B + Q = —4.0, G = 0.0, a = 0.3 and ~ = 0.2. As can be seen, the bias is large for events of smaller magnitude and also quite large for events around the detection threshold. In this case it is

30

35 40 45 True magnltude(m)

50

m

55

Fig. 2. Expected log(A/T) at a station for given true magnitude. Station parameters used are G 0.0, B + Q = —4.0, a = 0.3, 7 = 0.2.

241

not record the event were operating properly, the information that these stations did not record the event is useful and should be utilized in the estimating procedure. This was first done by Ringdal (1976), who considered the case with y 0, i.e., where the stations detect an event as soon as the “observed magnitude” rn1 is larger than G. Ringdal’s results differ from those given here, because he included in the sample space the cases where the event was not seen at any of the stations. This affects the estimated magnitude for smaller events, but for larger events and/or a larger number of stations, the differences will be small. The reason for this is that the probability of not seeing an event tends to zero as the true magnitude increases, and/or the number of stations increases. Suppose that we have a network ofM stations recording an event, and also that the event is seen by at least one station. The corresponding log likelihood is M

log L

=

~II~ log h,(a,/m) t—i

r—(B. + Q• + m

M

— log 1

—

FT I

=1

F

[

—

2

G~)1

j

(23)

‘~/cJ, + 7

We can see here that the likelihood is not a product of “individual likelihoods”. If the cases in which events are seen at all stations, or at a specific station had been considered, the likelihood (19) would have been a product of “individual likelihoods”. It is shown in Appendix 1 that the likelihood estliriator for magnitude is consistent and its asymptotic standard error is given. Can we use the conditional approach to estimate the parameters B,, G,, a, and 7, and magnitude simultaneously? To try to answer this question we assume that we have a sample of N events such that each event is seen and recorded by at least one station. In terms of “true” log(A/T) — Q arriving at the stations, we have the following model y,1—B,+rn1+e,1

i1,2,..,M;/1,2,...,N

(24)

If we had had the case where every event was always seen at every station, this would have been the standard model for the analysis of variance of a two-way classification. In order to identify the model we have to impose a normalizing condition on the B. (or m1) values. The reason for this is that these unknowns are identi-

fled only up to an additive constant. We will here adopt the condition ~ B. = 0, that is, the average station bias is set equal to zero. In the analysis of the variance model with fixed effects, the likelihood does not have a maximum if we treat the a~ values as unknown. This will also clearly be the case when we have missing observations, i.e. when some stations fail to report an event because of the detection properties of these stations. However, even if we treated the a, terms as known constants, there are difficulties because the number of unknowns increases with sample size. That is, as the number of events increases, the number of unknown event snagnitudes increases. Therefore, the standard methods for investigating the properties of the maximum likelihood estimator do not apply. The reason for this is that the estimator for m• is not sufficient independent ofB1, G, and ‘y,. And as we cannot obtain a simple relation expressing the estimator for rn1 in terms of data and the parameters B1, G• and 7~,we cannot use the approach used by Christoffersson (1970) in connection with the estimation of factor loadings in the case of missing observations. Therefore, at least for the present, we have to regard the conditional approach as a heuristic one for estimating the station parameters.

TABLE I Estimates of station biasB1, and station detection threshold

Station

NAO RES HFS UBO KBL FFC BLC ALE FBC CHG COL FCC YKC

B,

Q,

Gj

0 0.00 0.38 0.00 —0.16 0.09 —0.03 0.20 0.05 0.05 —0.39 —0.21 —0.15 —0.21

3.7 3:7 3.7 3.7 3.7 3.9 3.7 3.7 3.7 3.7 3.7 3.8 3.7 3.7

—0.1 +0.5 +0.1 +0.1 +0.2 +0.6 +0.9 +0.8 +0.8 +0.3 +0.4 +0.7 +0.7

242 TABLE II mb values Event

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

Station LAO

MBL

NAO

RES

HFS

UBO

KBL

FFC

BLC

ALE

3.52 5.88

—

—

—

—

—

—

—

—

—

6.08

—

—

5.85

5.91

—

—

5.99 4.75

5.70

—

5.92 4.69

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

5.18

4.66 4.05 4.76

4.97

4.92

5.09

4.98

5.20

4.92

—

—

—

—

—

—

4.99

4.88

5.03

4.86

—

—

—

—

4.67

4.65

—

—

5.00 4.25 4.56 4.57 4.4’)

—

—

—

—

—

—

—

—

—

4.00 4.26

3.80 3.70

4.03

—

—

—

—

4.33

—

—

3.82 5.10 3.73 4.92 4.00 433

—

—

—

3.88 335 4.13 4.06 3.80 3.71 3.72 3.71 3.95 3.55 333 3.76 3.88 3.78 4.56 4.05 4.83 3.90 4.12 3.38 4.44 4.09 4.34 3.81 5.83 4.47 4.15 4.37

4.45

4.58 4.05 4.09

—

—

3.90 4.57

4.00 4.16

5.28 4.25 5.21 4.28 4.74 4.62 4.50 3.88 4.41 4.72

—

—

—

—

—

3.66

—

—

—

—

—

—

4.00

4.02

—

—

—

—

—

—

—

—

—

—

—

—

—

4.50 4.30

4.42 3.90

—

4.35

—

—

—

4.38

4.08 4.01

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.23

4.05

—

4.23

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

5.44

—

—

4.72

—

—

—

—

—

4.49

4.69

4.98

—

5.12

—

—

—

—

—

4.30

434 4.77 4.47 4.45

4.89 4.22 4.58 4.02 4.02

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.90 4.50 4.45

430 4.38 4.03 3.79 5.73 4.25

4.90 4.72 4.37

4.48 4.01

4.10

4.47

4.87

4.48

—

—

—

4.47

4.42

—

—

—

—

—

—

—

5.63 4.07

5.58 4.18

5.96 4.41

6.07 4.50

5.86 4.50

—

—

—

—

4.43

6.06 4.50 4.57 4.64

4.05

4.44 3.90 5.68 4.10 4.07 4.20

4.89 4.46 4.25

—

—

—

—

—

—

—

—

—

—

—

—

3.95 4.33

3.83 4.13

—

—

—

—

—

—

—

4.30 4.69 4.08

—

—

4.35 4.62 4.40

—

—

—

—

—

—

—

—

—

—

—

—

—

—

5.19 4.47 4.77

—

4.30 5.83 4.59 4.90 4.28 4.20

—

—

4.86 4.01 4.55 —

—

4.97 4.37 4.61 4.00 4.26

—

—

3.90 3.92

—

4.66 4.26

—

6.17 4.62 4.38 4.47 4.20 4.20 4.57 4.40

3.41

—

—

—

—

—

—

—

—

—

—

—

4.08

4.36

—

3.90

4.01

—

—

—

4.55 3.78

3.88 3.89

—

—

—

—

—

—

3.93 3.45 4.04

3.95

—

—

—

—

—

—

—

—

—

—

—

4.09

4.50

4.68

—

— — —

—

4.65

—

3.80 4.18

—

4.68

—

3.77 4.20

—

4.10

3.92 3.55 3.87

243

Average mb FBC

CHG

COL

FCC

YKC

—

3.92

—

—

—

5.63

—

5.20

—

—

—

—

—

—

435 4.37

—

—

—

—

—

—

—

—

—

4.79

—

4.47

4.76

4.67

—

—

—

—

—

4.99

—

430

4.74

4.80

—

—

—

—

—

4.49

—

4.28

—

—

—

4.37 4.42

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.39

5.31

—

—

—

—

—

—

—

4.82

4.20

4.17

—

4.77

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.78

4.35

4.24

4.74

4.46

—

—

—

—

—

—

—

—

—

—

—

—

—

5.28

5.77

5.77

—

—

5.74

4.54 5.19

—

—

—

—

—

—

—

—

—

—

—

4.08

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.04

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.02

—

—

—

—

—

—

—

—

—

—

4.62

4.04 4.16

-

Max. likelihood mb

Standard error

3.72

3.55

5.80 4.66 4.37 3.82 4.93 4.01 4.92 4.27 4.51 4.31 4.26 3.71 4.04 4.26 3.73 3.91 3.72 4.21 4.11 3.55 3.53 4.07 3.88 4.04 5.00 4.28 4.68 4.17 4.22 3.38 4.57 4.36 4.35 4.01 5.75 4.36 4.29 4.32 4.20 4.13 4.35 4.29 3.41 3.88 4.13 3.78 3.95 3.72 4.33

5.25 4.06 3.60 3.24 4.87 3.67 4.86 3.94 4.43 4.11 4.01 3.25 3.90 4.08 3.39 3.70 3.12 3.94 3.86 2.83 2.79 3.81 3.30 3.56 4.67 3.98 4.64 3.81 3.93 2.45 437 4.08 4.14 3.87 5.76 4.22 3.90 4.16 3.32 3.83 4.16 3.73 232 3.20 3.92 2.65 3.76 3.71 4.25

0.157 0.078 0.101 0.148 0.244 0.079 0.138 0.079 0.019 0.085 0.097 0.104 0.240 0.112 0.099 0.192 0.133 0.306 0.109 0.116 0.619 0.698 0.120 0.221 0.154 0.081 0.106 0.081 0.121 0.109 2.354 0.082 0.099 0.096 0.115 0.078 0.092 0.112 0.095 0.214 0.118 0.095 0.129 1.763 0.260 0.111 1.009 0.127 0.131 0.091 (continued)

244 TABLE II (continued) Event

Station LAO

MBL

NAO

51

4.13

4.47

52

4.24

—

4.07 4.23

53 54

4.12 4.69

—

4.79

55 56 57

4.37 3.91 4.75

4.71 4.23 4.97

58 59

3.56 3.43

60 61

3.98 4.71

62 63

4.31 4.13

64

RES

HFS

UBO

KBL

FFC

BLC

ALE

—

3.83

—

—

—

—

—

4.20

—

4.25

—

—

—

—

4.41

4.60

4.59

3.95

4.48

4.28

4.58

4.39

4.89

3.93

4.34

4.47

4.54

4.75

4.50 4.66

4.77

4.71

4.79

4.38

4.93

4.62

4.84

4.70

—

4.59

3.93

—

4.29

—

4.50

—

434

5.02

4.71

4.24

4.55

4.64

4.80

4.61

4.08

—

4.12

—

4.02

—

—

—

—

—

—

—

—

—

—

—

—

—

4.45 5.02

4.08 455

—

—

—

—

—

—

—

5.04

4.19

4.46

4.43

4.63

4.80

4.66

—

—

4.57

—

4.27

—

—

—

4.22 5.42

4.29

—

—

—

—

5.14

5.58

5.42

65

3.78

—

66

—

—

67

3.50

3.44

68 69

4.03 4.28

4.29 4.40

70 71 72

4.04 4.32 4,67

4.28 4.05 4.68

5.15 3.60

5.56

4.61 5.13

4.41 4.82

4.60 5.51

—

—

—

—

—

—

—

3.88 335

4.15 4.08

3.70

—

—

—

—

—

—

—

—

—

—

—

—

—

—

3.92

—

—

—

—

3.91

4.76

—

—

—

—

—

—

3.78 3.32 4.36

4.42

3.75

—

—

—

—

—

—

—

—

—

—

—

—

—

4.45

4.33

4.73

4.49

4.98

4.50

For the case of estimating the magnitude only, the likelihood estimator seems to be fairly insensitive to moderate changes in the preset parameters. However, it is important that the stations that do not report an event have been operating properly at that time. Otherwise, there will be large biases in the estimated magnitudes. The asymptotic standard errors, on the other hand, are directly related to the present level of scattering, i.e., on the a~values and also to some extent on the 7~values. Tables I and II and Fig. 3 show an application for a simulated network consisting of 15 stations to an aftershock from Japan. Here all a, values are present to 0.3 and all 7~to 0.2. In the Figure there are some outliers. The first, event no. 2, illustrates the importance of the above-mentioned condition that the stations have to be in operation, and/or that no extremely high noise levels or interfering events are present. This is because this event is so large that it should have been seen by all the stations. Therefore it is probably best to have a variable threshold, for example related to the noise level, in the time interval

—

preceding the signal, as suggested by Ringdal (1976). The other outliers, events 3, 4, 40 and 47, are small events seen by 3, 1, 1 and 1 stations respectively, and reflect the bias in the conventional method for magnitude estimation. Average mb

6

2

3 6 ® 4~

3 .

4

~ .

.

.

hkelihood mb

Fig. 3. Average m~,and maximum likelihood mb for a simulated network of 15 stations. Aftershock consisting of 72 events.

245

FBC

CHG

COL

FCC

YKC

—

—

—

—

—

—

—

4.85

—

—

4.52 4.71 4.79

4.20 3.85 4.54

—

—

—

—

—

4.47 4.64

—

—

—

—

4.82

—

—

4.60

4.41

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.60

—

4.24

4.47

4.56

—

—

—

—

—

—

4.47

5.34

5.03

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

—

4.19

4.56

—

—

4.39

—

3.2. Unconditional maximum likelohood estimation Consider firstly the case with just one station. From eq. 14 we get, lettingrn0 -+—°°

G(a)da

/f

=

f

h(a/m)da exp(—j3m)dm

f

[1

—

Standard error

0.116 0.100 0.089 0.085 0.081 0.109 0.082 0.138 1.473 0.123 0.082 0.109 0.096 0.078 0.183 0.150 0.198 0.128 0.110 0.120 0.127 0.085

F(aN/m)] exp(—~3rn)dm

—

=

f 4[(B

2 +Q +m =

—

+72]

G)Wa (27)

const

Thus, the distribution of recorded events at the station [1 —F(a EN/rn)] exp(—j3m)d,n

(25)

is proportional to exp(—~3a) [(a G)/~] di

Now

(28)

—

Evaluating the proportionality constant (see Appendix 2) we fmd that the distribution is

h(a/rn) exp(—flm)dm

=

3.85 4.07 4.29 4.44 4.63 3.94 4.58 3.67 2.57 3.79 4.55 3.94 4.13 5.13 3.42 3.59 3.37 3.74 3.92 3.81 3.75 4.42

X exp(—j3m)dm

—=

f

4.12 4.35 4.38 450 4.67 4.24 4.67 3.94 3.43 4.17 4.60 4.24 4.26 5.21 3.69 3.91 3.64 4.08 4.34 4.05 3.97 4.56

4.57

5.39

—

Max. likelihood rn,,

—

3.81 353 5.00

—

Average

f

G(a)da = exp(J3G 2/2a2}

(1/~~J~ a) exp{—[a

X exp(—~rn)d~’[(a X exp(—j3a) 4 [(a

—

—

—

(B +

G)/71 drn G)/7]

=

Q

+ m)J

const (26)

—

~y2j32/2) exp(—g3a)

4 [(a

—

G)/

7] (29) It should be noted that thi~distribution is independent ofthe scattermg standard deviation, the region— station bias, and distance—depth correction, i.e., the observed distribution depends only on the seismicity parameter 13, and the station detection parameters G and ~.

246

Turning to the likelihood estimator, we have for a sample of N independent observations, a1, a2, aN

Although not explicitly stated, these likelihtod equations (3 1—33) were obtained by Kelly and Lacoss (1969), who in addition were estimating the total number of earthquakes in a given time period., The likelihood estimator for /3 considered by Aki (1965) can be obtained as a special case from the above equationsbyputting7=OandG=min(a1,...,a~). To ifiustrate the above method data in the distance range 30°—60°observed at the stations UBO, TFO and WMO The (see results Figs. 7—9 for locations) been analyzed. in terms of base 10have logarithms are shown in Figs. 4, 5 and 6. These show that the model fits the data reasonably well. However, for UBO (Fig. 4) there is a statistically significant deviation between model and data. On the other hand, the number of observations is large (4562), and hence so is the probabifity of detecting small deviations. The maximum deviation between data and model is 0.02. Whether this deviation is of any practical significance

...,

logL = Nlog /3 +N(IG N

—

i 2 ~N7 (32

N

a. + 11

E

/

.

—

log 4(’~” G)

(30)

7

‘~

1=1

The derivatives of log L with respect to /3, G, and 7 are 2/3

—

N a,

~ log L/~j3= (N/j3) + NG

—

~ log L/~G= N/3

~[4’y,)/cb(y,)] i=i

(31)

11

N7

N —

(1/7)

(32)

N

~ log L/&7 = Nj32’y

—

(1/7) ~ y,[d’(y,)/F(y,)]

(33)

—1

where y,

= (a, — G) 7. The likelihood estimator is defined as the solution

wifi depend on the situation in which the model is applied. For example, in seismic risk studies this

to

deviation may be of very great importance. Considering the estimated /3-values we see a rather large variation, ~UBO = 0.85, ~WMO = 1.15 and ~TFO = 0.93. This difference may be attributed to the differ-

log L/613 = ~ log L/~G= ~ log L/S~= 0 which corresponds to the maximum of log L. t~

(34)

Frequency

1.0’

0.9 0.8 0.7-

significant deviation (a

.05)

0.60.50.40.3-

0.20.1-

0.0 I

I

I

I

I

I

I

3.2

3.4

‘—1og 10(7IJ5’)

0.0

0.2

0.4

06

0.8

1.0

1.2

14

tO

1.8

2.0

2.2

3.4

2.6

2.8

3.0

3.6

Fig. 4. Station UBO. Distance 30—90°.4562 events. Observed and theoretical log jo(A/T) cumulative distributions (3 = 085(0.017), a= 0.19(0.010), 7=0.11(0.007).

247 Frequency 1-0. 0.9• 0.8-

0.70.605 0.4 0.3

0.20-I,

o~o’

I

-2

I

I

I

I

I

I

I

I

I

I

I

I

0.0

0.2

04

0.8

0.8

1.0

12

14

16

13

2.0

2.2

I

24

I

I

I

2.6

2.8

3.0

Fig. 5. Station TFO. Distance 30—90°.3725 events. Observed and theoretical log

1o(A/T) cumulative distributions /3 = 0.93(0.019),

0=0.09(0.008), 7=0.10(0.005).

ent seismic areas (see Figs. 7, 8 and 9) that are sampled. The southern part of Japan is outside the distance range 30°—90° for WMO, which may explain the rela-. tively high a-value for this station. On the other hand,

it has been suggested that the linear seismicity model is valid only for a limited magnitude range, and that the slope of the seismicity curve increases with magnitude. The relatively high detection threshold and

Frequency

1.0’

0.9 0.8

0.7

0.6 05

0.4 0.302

0.10.0—

I

0.0

0.2

0.4 0.6

I

I

I

I

I

1

0.8

1.0

12

14

16

1.8

I 2.0

I

I

I

I

2.2

24

2.6

28

33

iog

10

(I~/’F)

Fig. 6. Station WMO. Distance 30—90°. 2321 events. Observed and theoretical log10(.4/T) cumulative distributions (3 0=0.40(0.015), 7=0.16(0.007).

1.16(0.035),

248

UBU UINTA BASIN RR.

7

Fig. 7. Azimuthal projection with center at station UBO. The distance curves from origin are 30,45, 60, 75 and 90°. Azimuths are

0, 30,60 and 90 (and soon) degrees.

/3-value for WMO support this later hypothesis. Whereas the magnitude distribution’for one station can be used to obtain reliable estimates of the detection parameters G and ~, no information about the station bias B and the scattering a can be obtained, as

the distribution of observed magnitude at a station is

independent of these parameters. To estimate these parameters we would like to consider the distribution (14) for the general case of M stations. Unfortunately, evaluation of the integrals in (14) for the general case

249

WMO

WIC~-~IT~ MTS.

!~—;~

7

/

-4ç~

~

‘

/

~ ~‘

~

~

‘~%~ ~

~ ~

—*~

~

k

-° ‘S.

I~.

-

.

1~

-

a ~

‘I Fig. 8. Azimuthal projection with center at station WMO. The distance curves from origin are 30, 45, 60, 75 and 900. Azimuths are 0, 30, 60 and 90 a.s.o. degrees.

leads to very complicated expressions which do not seem to be useful from a practical point of view (except

integration does not seem to work with more than a few stations.

of course for the case in which events are seen jointly by all stations, but then it is very difficult to obtain enough events for analysis). The use of numerical

However, for the case M = 2 and events seen jointly by the two stations, the unconditional distribution is not too complicated. In this case the distribution (see

250

TFO

TONTO FOREST AR. ‘I

~k

~

7

~ \ ~

1

1

/

c...

Fig. 9. Azimuthal projection with center at station TFO. The distance curves from origin are 30, 45, 60, 75 and 90°.Azimuths are

30, 60 and 90 a.s.o. degrees.

Appendix 2) turns out to be

where

G(a1a2)=(1/V~/~) G 1)/’y1] cli [(a2 G2)/72] exp(—y/2)

y

2

[a1—a2 —(Bi ~Qi —B2

—

Q2)]

—

+/3[2a~(ai —B1

—

I

[4(ai) exp(a2) + 43(a~)exp(a4)]

(35)

—

—

Qi)+2o~(a2—B2

aWfl)]/(a~+ o~)

—

Q2)

(36)

251

[(B1 + Q1

=

/

a

Sv/0~ +

B2

—

Q2)

—

—

(G1

—

G2)

—

(c4 + ‘y~)/3]

c4 + 7~+ 7~

(37) (38)

2(a~+’y~)/2

following distribution g(a2/a1)da2 2}

1

exp~—[a2—a1

2—j3(B2 +Q2—G2)+j3

—Q2)—B]/2a

—(Q~

[(B

0.3 =

2 + Q2

B1

—

Qi)

—

—

(G2

—

G1)

—

(o~+

~ ~,(a2

+ ~~‘+‘y~ 2(al +7~)/2 /3(B1 +Qi G1)+/3

0.4

(39) (40)

—

When comparing two stations it is a common practice to plot the magnitude at one station for given magnitude at the other. It is shown in Appendix 2 that the conditional distribution of log(A/1) at station 2 for given log(A/2’) at station 1 is g(a2/a 1) = 4) [(a2 G2)/72 1(1 /‘s/~ a) —

X exp [a~ a1 4){[ai +B + {—

/

—

—

—

(Q2

(Q 2

where B=B

—

Q1)—

Qi)

B] 2/202 2+ 7~} (41) G2]/’~/a —

}

(42)

2i

2 B1 a2 = o~+ a22 —

—

/30

(43)

2/./~

E(a2/ai)

=

B + (Q2

—

+ a1 +

Q1)

=

7~)

X (4)’(x)/4)(x))

(44)

and D2(a

4/(o2

+ 7~)](4)’(x)/4)(x))

2/a1) = a2 [o X [x + (4)’(x)/(4)(x))]

(45)

where B and a2 are given by (42), (43) and x = [a 2+ (46) 1+ B + G2]/s/a If we had also considered events not seen at station 2 (but seen at station 1), we would have obtained the —

Q1)

—

Qi)

—

—

G2] /~/~‘~

(47b)

}

if not seen at station 2.B anda are given by (42) and (43). This distribution is of the same type as (6). Ringdal (1975) based a likelihood estimator for the detection parameters on the probabilities of detecting an event at station 2 for a given observed magnitude at station 1, i.e., essentially on eqs. 47. However, with this latter approach it 2) is not the and possible slope ofto theseparate detection effects(72). of scattering (a the method based on distribucurve To illustrate tion (41), the same data as for the single-station case have been used. The analysis has been made for both log(A/7) and magnitude data. Letting * denote the parameters in the model with magnitude data, the relation to parameters with log(A/7’) data is B

=

G

0 2 + ~2; a

=

a;

=

~

+

S~ 2

(48)

where ~ is the average distance—depth correction and S~2is the corresponding variance. The results of the analyses are shown in Tables 3 and 4. In Fig. 10 a plot of the observed average magnitude at WMO for given magnitude at UBOmagnitude is shown, calculated together with expected average fromtheeq. 44. There were no difficulties in solving the likelihood equations except for one case where log(A/7) data were used to estimate the parameters for TFO, with WMO as the reference station. This problem was probably caused byof insufficient coverage of events in the neighborhood the detection threshold for TFO,

—

(Q2

+ B(Q2

~{

0 = B; G

(a

(47a)

/

if seen at station 2 and g(a 2/a1)da2

The parameters that are identified are B, G 2, a2 and ‘y~,and can thus be estimated. The corresponding likelihood estimator is given in Appendix 2. l’his distribution is of the same type as the distribution of observed log(A/7) for given true magnitude (eq. 20), so the first two moments of (41) are

G2\

7

______

~

—

due to the high detection threshold for the reference station. Regardless of this problem, by restricting the permissible range of variation for the estimated tineshold so that it had to be larger than the smallest observed log(A/7) + 0.01, it was possible to obtain reasonable estimates of the parameters. Comparing the results based on log(A/7) and magnitude, we find close

252 In the distribution (41), it is not possible to iden-

TABLE III Log (A/i’) data. Standard errors given in parentheses Station

Reference station

B

G

a

WMO

—0.04 (0.01)

0.07 (0.01)

0.37 (0.01)

0.05 (0.01)

0.15 (0.01) —0.34 (0.02)

0.11 (0.02) 0.37 (0.02)

0.37 (0.01) 0.37 (0.01)

0.08 (0.02) 0.17 (0.01)

—0.07 (0.02)

0.39 (0.02)

0.35 (0.01)

0.18 (0.01)

—0.44 (0.01)

0.03 (0.01)

0.37 (0.01)

0.08 (0.01)

UBO TFO UBO

tify either the individual scattering variances or the region—station biases. However, if it can be assumed that the parameters are at least approximately constant over the different regions, the results for all possible combinations of station pairs can be combined to obtain estimates of scattering variance and region— station bias. Using log(A/7) data we have for the different combinations (see eq. 43) 4Fo

=

0.372, 03~TMO+ 4F0

=

+ U~rMO= 0.372, 0~JBO+ ~

+

=

0.372 0.352

WMO TFO UBO TFO WMO

—0.29 0.01 * 0.36 0.05 (0.01) (0.02) (0.01) (0.02) ___________________________________________ * Lower limit for G.

~~JBO + aTFo =

0.372,

0~%TMO+ OIFO =

0.362

Solving this system by least-squares, we obtain the following estimates of the scattering standard deviations: = 0.27; ~WMO 0.25 and UTFO = 0.25. Similarly, for magnitude data we get UUBO = 0.27; ~WMO 0.25 and UTFO = 0.24.

-

Turning to the region—station bias, we obtain from

eq. 42, using log(A/7’) data agreement. If we used eq. 48 and transformed log(A/T)

Bu 00

parameters to magnitude or vice versa, the results would differ only slightly. There is a tendency for magnitude data to give somewhat larger numerical values for both G and

—

BWMO —

Bu~o BTFO —

—

I3OIJBO =

/3G~JBO=

—0.04,

0.15

~.

TABLE IV Magnitude data (standard errors within parentheses) Station

Reference

(49)

B*

G*

WMO

—0.05 (0.01)

3.82 (0.03)

0.37 (0.01)

0.21 (0.02)

3.72

0.17

TFO

0.14 (0.01)

3.88 (0.03)

0.36 (0.01)

0.24 (0.02)

3.71

0.16

UBO

—0.32 (0.02)

4.17 (0.04)

0.37 (0.01)

0.24 (0.01)

3.75

0.17

TFO

—0.06 (0.02)

4.20 (0.04)

0.35 (0.01)

0.26 (0.02)

3.74

0.16

UBO

—0.42 (0.01)

3.77 (0.03)

0.36 (0.01)

0.22 (0.01)

3.71

0.17

WMO

—0.26 (0.01)

3.67 (0.03)

0.35 (0.01)

0.20 (0.03)

3.71

0.15

a*

y*

SQ

station

UBO

WMO

TFO

253

ing to BUBO and BWMO. Using magnitude data, we obtain EUBO = 0.28; EWMO = 0.12 and = 2.52 (in

a

1 .8

base 10 logarithms, ~ = 1.10). The reliability of this a-value is of the same order as for log(A/7) data. It must be noted that the results concerning bias and scattering are based on the assumption that the parameters are the same for all sampled regions, and that analysis of single-station data indicated that the slopes (J3) may not be equal. This difficulty can be overcome by using data from one specific region. However, in this paper, data are used only to illustrate the statistical methods and are considered to be of suffi-

/

/

4

:~ 60

~

:~

/

/

‘6

X

X

x x ,~

‘7

cient quality.

4

4. Conclusions 4./~

321 6i4b1 4

81 21

e eo’

.~I ~

12

.8 2 6 5.0 .4 .8 2 .6 70 4 8 Fig. 10. Average observed magnitude at station WMO for given observed magnitude at station UBO for 1958 events in the distance range 30—90°./3~—0.32(0.02), a 4.17(0.04), S 0.37(0.01), 70.24(0.01), Q 3.75, SQ 0.17. 4

=

BWMO

—

BUBO

—

(3G~VMO= —0.34,

BWMO

—

BTFO

—

I3°~wMO = —0.07

BTFO

—

BUBO

BTFO

—

BWMO

—

—

=

=

/3~iFO= —0.44, o~FO= —0.29

The conventional estimator for network magnitude is an average of the observed magnitudes of the detecting stations. This may result in seriously biased estimates, especially for events with magnitudes close to or below the 50% detection threshold. By using the conditional maximum likelihood estimator given here (or the one given by Ringdal (1976)), this bias can be effectively removed. Joint maximum likelihood estimation of the detection, region—station bias and scattering parameters for a network is very complicated. However, by considering all combinations of pairs of stations in a network it is possible to obtain reliable estimates of these parameters.

Appendix 1 (50)

The conditional distribution of observed log(A/7) at a station is given by eq. 19 as

We first note that the region—station coefficients can

only be determined up to an unknown additive constant, and therefore adopt the normalizing condition BTFO = 0. If in addition we substitute the estimated scattering variances obtained from (49) for 0?JBO, a~4()and 4F0~we can solve (SO) for BUBO, BWMO and (3 by least-squares. This gives EUBO = 0.31; BWMO = 0.14 and = 2.78 (in base 10 logarithms, = 1.20). This estimate of /3 must be regarded as very unreliable compared to the estimates ofBuBo and BWMO. This is because the diagonal element in the inverse of the moment matrix obtained when solving (50) is more than 100 times the elements correspond-

a

a

H(a/m)

=

X exp

(1/v’~a)

{—

[a — (B +

G + m)} 2/202 }4)(y)/4)(x) Al 1 ______

wherey = (a G)/’y,x = (B + Q + m G)/~/o2+ ~y2, and 4)(x) = (1/v’s) exp(—t2/2)dt. To obtain the expectation of observed magnitude for given true magnitude we first evaluate —

—

f~i..

F

=__~

~

J’ —~

[a

—

(B + Q + m)]

254

X exp{—[a

—

(B + Q +m)]2/2a2}4)(a _G)~ (Al .2)

2+72~magnitude for To obtain variance of observed wherex (B the + Q +m G)/~/a —

given true magnitude we first evaluate

Weputua—(B+Q+m)andG 0—B+Q+m—G, thus

F 2/2o2)4)[(u + G

f

F (l/~’~a) u exp(—u

0)/7] du (A1.3)

2/2a~)= (—u/a2) exp(—u2/2a2) d ~ exp(—u and so

f (l/~./~a)(l/a2)[(u2/a2)

1] 2/2a2)4)[(u ÷ G X exp(—u 0)/-y]du 2/2a2)4)[(u + G = {(-u/~/~a~) exp(-u 0)/7] -

(Al .4) +

f

(u/2ira3~’)exp(—u2/2a2)

—=

F = f(—a/’.~/27r)exp(—u2/2a2)4)[(u + G

2/2’y2]du = f(u/21Ta37) 0)/y]}°~

f

+

(a/-~J~i) exp(—u2/2a2)(l/.J~7)

0)

+ (u + Go)2a2]/2a2’y2} du

—G

2 + 72)]/[./~~,/o2 0 exp[—G~/2(a [—x/(a2+ ~2)]$’(x) withx as before. This gives =

o+f(a/27r7)exp{_[u272 (A1.5)

Now

E(u2/m) =

a2

—

[a4x4)’(x)/(a2 + 72)4)(x)]

+ 72(02 + 72)]

(A1.ll)

(Al.l2)

and

u2”? + (u + G 2a2 0) X [u +G 2/(a2 0a We then find

D2(u/m) = E(u2/m)

= (02 + 72)

2 +~2]2

exp[—G~/(a2 + 72)2]

+(G~/o2 2) 7

(Al.6)

J’ 2~ry

D2(a/m) =

~

0o

(Al 7)

and finally 72

E2(u/m)

(a~+72)4)(x) __________

I

4)(x))

(AI.13)

+____

~2

a4

—

4)’(x)

(a2 +72) 4)(x) It can be easily shown that

)

/ lim 4)(x) (x \ + —

=

(

x+—) 4)’(x)’ (A1.l4) cli(x),

(Al .15)

0

and

______

(a2/../02 +

U

—

a4cli’(x)

Thus, the variance in the observed distribution is

X exp{—(a2 + 72)[u + G 2/(a2 + 72)12/20272) du

F

X exp[—(u + G0) X exp(—u2/2a2) exp[—(u + G 2/2’y2]du

X exp[—(u + Go)2/2’y2]du

F

~

\/3J~)

exp[—G~/2(a2 + 72)] (Al.8) lim

Then E(u/m) = F/4)(x)

[4)’(x)/cli(x)]{x +

[4)1()/4)~]}

=

(Al .9)

and

(Al.16)

This gives lj.m D2(a/m).a2

(Al.17)

m-~

E(a/m)B+Q+m+F/4)(x) B + Q + m + (a2/~/~?+ 72)[4)’(x)/4)(x)]

1

X~O~

and (A1.l0)

lim D2(a/m) =

m -~-

~2

[1 — 1/(1 + 72/02)1

(Al.l8)

255

Turning to the likeithood estimator we have

asymptotic results for maximum likeithood do not apply directly. However, by showing that the maximum for L(m) is the same as for L (m) when M tends

log L = log H(a/M) =

const

[a

—

—

(B + Q + m)]2/2a2

—

to infinity, weiscan apply likelihood. the asymptotic formulae to L0(m) which a regular

log cli(x) (A1.19)

8 log L/8m = [a —

—

We have

(B + Q + m)] /a2

log L(m) = log L 0(m)

(1 /s,/02 + 72) [4)’(x)/4)(x)]

(Al .20)

a = B + Q + m + (a2/V’a2 + 72)[4)’(x)/~x)] (Al.2l)

M

L(m)

=

‘

M

H hi(ai/m)/(l

0(m)/P0(m) =

~ (Al .23)

L

=

1

(Al.3l)

=

M-~ It then follows directly that

Jim log P0(m) = 0 M

~1h,(a,/m)

(A1.24)

M —

Jim

8ii

M—=

£1

P0(m) = 1

[1 F(a,eN,/m)

(Al.25)

1=1

1

=

(Al .32)

/

2a1 }4) [(a,

log P0(m)

Oforallfinitemandn= 1,2,3,...

8m

(Al.33) From this it follows that the asymptotic properties of the estimator obtained by maximizing L(m)are the

exp{—[a,

+ B- +

—

—

(B, + Q, + m)] 2

G,)/y~]if the event is seen —

Gi)) if the event is not seen

same as those obtained by maximizing L (m). As the latter is a regular maximum likelihood estimator we have, letting m 0 denote the true magnitude of the event and m the likelihood estimator, that m0) tends to a normally distributed variable with expectation zero and variance 0(m) 2 -1 D2[~mmo)]([~E( alogL 8m )/mmo]}. ~M(m

—

(Al .26)

(Al .34) Now

and _______________ P(a,eN,/m)

(Al .29) (A1.30)

M—*

L0(m)-~

h,(a,/m)

1 —(1 _&)M~
—11 P(a 1eTV1/m))

with

(Al .28)

L(m) for every finite m. Specifically, we assume that there exists a ~> 0 such that (a,+B 1+Q1—G,)ij~l_6 ~0~+77 foralli

_____

2Wa2+72)z(x) (A1.22) a-B+Q+m+(a Turning to the case of several stations, the likelihood is (from eq. 19)

log P0(m)

Provided B,, Q,, G,, 01,71 are well behaved for all i, that is, are such that each station has a non-zero probability to see the event, it follows that log L(m) tends to log

Solving 8 log L/8m = 0 then leads to solving

or, letting z(x) =

—

Im

=

m

=

4)(_(m +B, ++Q~ Gi)) —

(Al.27)

It was mentioned earlier likelihoods, that this likelihood is not a pro~uct of independent so the usual

(l/M)Et(8 lo~L*(m)yI 8m I } 2/m = mo) M = M 1 ~,=i EI(8 log h,(a,/m))

(Al .35)

256 a

1

~B, + Q~+ m0)

—

Appendix 2

a?

The unconditional distribution of observed log (A/fl at a station is (from eq. 23) proportional to 0(a) = exp(—13a)4)[(a G)/ G 7] (A2.l)

if the event is seen

—

—l log h,(a,/m)’ I

As the integral of the distribution equals 1, we can

\/‘O? + 7?

r(0

(rn0

L

8m

)/m

=

rnol

x[4)’(

=

+B, +

—

G,)\ determine the proportionality constant by evaluating

—

/ 4)l(

(rn0+B~+Q, ~o?+7/

G~)) ~)i

-

I -=

exp(-/3a)4)[(a

G)/7] da

-

X exp(—(3a)4)[(a

G)/7]

—

=

}~00

{(-1//3) +

if the event is not seen (A.136)

x f exp(-~)(l/~7)exp[-(a

0)

(Al.37)

1 f =o+~

G,)/-%/a? +

(Al.38)

2/2y21 da -

—=

or letting

G)

u,a,—(B,+Q,+m x,

=

(m 0 + B, +

Q, +

we have

2/2 2} 7

V~7 exp{—[a—(G—2/37)]

X exp[—/3G + (3272/2] da

R

=

(1/fl) exp(—(3G + /3272/2)

(A2.2)

(ui/a?

)/~

1/8 log h,(a,/m)

=

mo] =

—

I

4)’(x,)

Thus, the distribution of observed log(A/7) is

I. ~sJü/+ 7~4)(—x,) ______

G(a) = j3 exp(f3G

/3272/2) exp(—/3a)4)[(a

—

—

G)/y] (A2.3)

(Al.39) and

E

2!

1

=

(a log 8m hi(ai/m)) 1mm I =

f

Suppose we have a sample of N independent events parameters 7 is obtained by maximizing the of the recorded at /3,theG,station. The likelihood estimator

logarithm of the likelihood. We have, letting a denote the observed log(A/fl

(1/~/~aXu?/a4) exp(-u2/2a?)

4)(Ui +

G 0\ Jdu,

“

/

N logL(/3, G, y)~I~logG(a,)=Nlog(3+N(3G i=1

+ [l/(a? + 7?)] [4)’(x,)/4)(—x,)] 2 4)(—x,) (Al .40)

2/32

From eq. Al .11 it follows that E({ [8 log h,(a,/m)/8m] /m = rn 2) 0)

I 1 \ r4)’(X,)11r4)’(~,)_1 2 \01/ Ic,

÷

— ~

1

—

~II~ 11

/3 N a• + N log ‘li[(a, 1=1

Ny

G)/7] (A2.4)

The first-order derivatives are, lettingy,

7?itL4)(—x,)i_x13 =b,

=

(a,

—

G)/7

N

(Al.41)

and finally

~

(A2.5)

8/3

/3

11

M

2 [.J/~i(m m —

D

—

0)]

=

~I b,)

l/(~

(Al .42)

8 log L(j3, G, ~ 8G

=

N(3

— -~-~

7 i=i

4)(yj) 4)~~j’,)

(A2.6)

257 alog L(/3, G, 7)

=

—N/327

87

— ~

N

/m+B

711 —

82 log L(J3, G, 7)/8/307 82 log L(J3, G,

=

N ~)=_~~

2

/1 ~

(A2.7)

Turning to the second-order derivatives we have 2) Ny2 (A2.8) 82 log L(j3, G, 7)/8(32 = —(N//I 82 log L(J3, G, 7)/8f38G = N (A2.9) —2N/l’y 4)’(Yi)(

X exp(—j3rn) dm

h,(a,/m)

1

=

1

8G8y

72

+

(A2.17)

Gi)

7i —

—

—

~(Y~)

—

Q,

—

—

Q1,

—

—

N •1

4)’(y,) (1 4)(~

~‘(yj)\

2

=

F1

2 log 872 L(j3, G,

~

=

—N/I

f

h1~a1/m)h2~a2/rn) exp(—/Im) dm

—=

f

b\ 4)(~)4)(~) 2 1 exp{—[(u ~ 1T01c72 + (v m)2/2a~] exp(—flrn) dm

=

+..LE N y,4)’(y,) 72 11

=

— ~1

—

(A2.12)

—

~? y, ‘li’(y,)/4)(y,))

(A2. 13)

—

2logL(0) 802

)ooo~j

1

=

2

=

f —=

L(0)1

Jo—00

h1(a1/m)h2(a2/m) exp(—/Im) dm

exp{-[(u

4)(rn+Bl

v)2

-

71 7~ \/\/a~+022 + /3(2a~u+ 2a~u— a~a~13)]/2(a~ + a~)}

F

f

(A2.19)

F

0 = 0] computed from the data instead of E[(82 log L(0)/802)/0 0~]because

Next, consider the case of events seen jointly by two stations. The joint distribution of log(A/T) values, a1 and a2 is

—

and we find

Next

(A2. 15)

m)2/a~] !3m

—

(A2.14)

802

—

+ f3(2a~u+ 2a~v a~o~j3)J

2 is given by eqs. or A2.8—A.2.13. where 82 samples log L/80 we may use (l/JV)[(82 log L(0)/802)/ For large

[82 log

m)2/a~+ (v

—(l/2)[(a~ + a~)/o~a~] (rn —+ a~)] ~[l/(a~ + — 2 [l/2(a~ [(u —+ a~)](o~u v)2 + a~v—a~a~j3)})

1

1 182 log L(0)] plim—I 2 N—’-°~N L 80 0—0~

(A2.l8)

=

—

Oo))= (_WE[~

rn)2/2o~

Now —

-

—

} —=

4)(,)

Let 0 = (j3, G, y) and let 0 = (J3, G, y) denote the likelihood estimator and O~= (/3~,G 0, 7o) the true parameter As wethat herethe have a regular likelihoodof estimator itset. follows limiting distribution O~)is normal with expectation zero and covariance matrix

G(a1a2) =

—

—

82 log L~,G, 7)

—

(A2.l6)

Xexp[—(a, B, rn)2/2o~] for i = a, 2 Let u a u a2 B2 Q2,a = a1 G1 and b = a1 B1 2 G2, then

4)’(Yi))

72

)

~~‘1th

(A2.ll)

X ‘(2

)

(A2.lO)

8G

02

1+Q1—G1\ (rn+B2+Q2_G2\ ~ + 7~ ~ ~ + 7~

(A2.20)

+Q1 _Gl)4)(m+B2 +Q2

—

G2)

2 -../a~+

~Ia 2+ (A2.21)

X exp(—f3rn) dm

Let —Q1—B1+G1_=c,—Q2—B2+G2=d,s1— + ‘y~,s2 = ~/a~ + ‘y~and s = ~/a~ + a~+ 7~+ 7~. Then

F2

=

f 4)(~_—_c~4)(rn_— \ S1

,

\

S2 d) exp(—(3rn) dm

258

1 ~4)(m_—_c’~ 4)(m_— d) exp(_/3m)1 13L \ S1 / +

1

m—c\ J ~—) =

The marginal distribution of a1, given that the event is seen at both stations is

\ S2

2/2s~]

I

~

exp[—(m

g(aj)f

d)

—

X exp(—f3m) dm

=

—=

J

G(a1a2)da2 4)(a2

—

exp(—{[a1 —a2 72

(

1

1

4)m—d’

+ (3~

)

~

X exp(—jlm) drn 2/2s~ jIm —

—[rn

=

(d

—

C)2I2S~1

-

52(3)12/2S2

—

—(rn d) —d13 + (l/2)f32s~

+

(A2.22)

-

Q1 -82- Q2)1

2

I 2(a~+ a~))da2 =

(.8k

2 -2(3a~(a2 B

—=

exp[-(rn

~

—

/

-

Q2)}

2a~(a

C0 exp~—[j3~a~2(3

1 B1

______

—

—

—

Q1)]

—

(A2.23)

/2(a~+ a~)}~

+

a~4)[(B 2 + Q2

Similarly —(m —

—

2/24

c)2/2s~—f3m c/I + (1/2)13252

—[m

—

—

(c

s~f3)]

—

(A2.24)

So

—

—

(3o~ G2)/~/~ + a~+ 722]

ai

—

(A2.32)

From this we can obtain the conditional distribution of a

~

2s~/2)4) F2

B1

+

=

exp(—j3d + f3

2, given a1 for events seen at both stations. We get —

S 2

—

C)

i~IC

—

s1

—

d\

g(a2/ai) = G(a1a2)/g(a1)

[4)(a2 —_G2’)

=

~ (A2.25)

X exp(—/3c + /32s~/2)

Xexp[—(a2

a1

—

—

1

/,,~

72

B)2/2a2]]/4)[(ai

+8

—

and finally G(aia2) =

_

G2)/~J~ +y~] (A2.33)

1

~

\ 4)(a1_—

\ 7i

+

Gi’~~4)(a2 G2~ / \ 7~ / —

X exp(—y/2)/[4)(a1) exp(a2)+ 4)(a3)exp(a4)] (A2.26) where 2 y = {[a1 a2 —(81 + B2 Q2)] + /I[2a~(ai B 1 + 2a~(a2 B2 Q2) —

—

whereBB2 + Q2 —B1 j3a~and ~2 = + This distribution is of the same type as that of ob—

—

served log(A/fl for a given true magnitude (eq. Al.l), so the first two moments ofg(a2/a1) are given by eqs. A1.10 and Al.l4. On the other hand, had we also considered the events were the events recorded notatseen station at station 1, it is2,readily given that seen

—

—

—

a1

a2 a 3

Q1)

—

—

a~o~f3]}/(a~ + a~)

=

[(B1 +Q~—B2— Q2)—(G1

—

(A2.28)

2(a~+ y~)/2 I3(B2 + Q2 G2) + j3 = [(B2 ~Q2 —B1 Q1)—(G2 G1) (o~+ y/3]/~,/~+a~+ + y~)

(A2.29)

—

—

—

a

2(a2 + y~)/2 =

13(Bi + Qi

—

I

G2)

-(a~ +y~]Wa~+a~+y~+y~ =

that we obtain the following distribution

(A2.27)

—

4

—

G1) + f3

I

g(a2/a1) da2

4)(a2

—

G2)

7~ -

1

exp(—{a2

2/2a2)da

[a~ +B])

2

=1 if seen at station 2 2 +~ L ~‘a

I

(A2.30)

if not seen at station 2

(A2.34)

(A2.31)

The log likelihood based on (A2.33) for a sample of N

)

259 independent observations is, except for a constant

2 log L/802/0 = O)J_1. Turning to these secondX(0 order derivatives, we find after some calculation that

N

logL~ —Nlogy— [a

2logL 21—Q21—~a11—Q1~)

O

N

0B2

a2

1 z~x,)Ex,+ z~x,)]

02+72j=1

(A2.43)

—B12/2a2 +log 4)[(a 21— G2)/72]

—log 4)[(a1, +B+ Q21 whereB(B Putting G = 2G

—

—

2 +y~] (A2.35)

G2)/sJa

2 = a~+ a2. 2 B~—$a~)and a 2 and 7= 72, Q~= Q2~ 2 +~2 (A2.36) Q~+a11 G)/Va

1 ~$+ y~= Ia

logL = ~N 8G2 1=1 +

—

—

x

02

Q1,

Z(y~)[yj +Z(yj)]/72

1 2 + ~2 z(x,)[x + z(x,)]

G)/y

—

(A2.37)

and z(x) = [4)’(x)/4)(x)]

(A2.38)

The first-order derivatives are 0logL_’~ ~a~ OB i=1 1 B ______

—

—

2 —

Q,

—

8a2

—

1

z~x,)

a

2

—

Q,

—

a21)

~=i + 72)2](a2 fx,[x, + z(x,)]

—

2} +

+ [x,z(x,)/(a 82

— ~

2 . 3/~4 21 B

—

2~

(A2.44)

a 82logLN~~ —~a

~2)

(A2.45)

log L = ~ yjz(y,) {2

0y2 ~ + Ex,z~x

7 2+

—

72)2(y2

+ z(y,)]

{x,[x, + z(x,)]

}

2} +

—

02)

1)/(a

a11)/a

(A2 .46) (A2 .39) N

0 log L OG

N =

~

—

z(x,)

1

z(y,)/y +

(A2.40)

02logL_~ z(x, OBOG a2 + ~2 [x, + z(x,)] —

(A2.47)

—

i—i N a2logLE —(a

2/a3

8 log 8a L

=

1=1 (az,

0

~

—

B

—

Q,

—

OBOa

a11)

—

—

y,z(y,)/y + z~x,)

(A2.42)

The likelihood estimator is then obtained by solving

=

=

(8 log L/OG) = (8 log L/8o)

Let 0~= (B, G, a, 7) denote the true parameter set and = ~‘) denote the likelihood estimate. It than follows directly from the general properties of the maximum likelihood method that i.ji~i(Oo is asymptotically normally distributed with expectation zero gnd covariance matrix given by [—(1iN) E(82 log L/802/0 0)]~which may be estimated by [—(1/N) (E,

O,

82 logL OBOy

yz~x,)

—

j=~—

(a2 +

~2)3/2

{x,[x, + z(x,)]

l} (A2.49)

—

N 82logL~ 0 {x,[x,+z(x,)] —l} OG8a ,=~ (a~+ 72)3/2

(0 log L/O’y) = 0

O

1)

—

(A2.48)

(~2+72)

(8 log L/813)

[az(x,)/(o2 + 72)3/2] {x, [x, + z(x,)]

(A2.41)

+ [ax,/(a2 + y2)]z(x,)

8 log L =

21—B— Qj—a11) 2/os

t=i

&,

— ~)

(A2.50) N

82 log L = 8G0y j=1

—

Z(y~)

+ z(y,)1

~

-y + [7z(x,)/(o2 + y2)3/2

]

—

{x. [x, + z(x,)]

1

}

—

1

}

(A2.51)

260

and finally 02_logL~ OaOy 1=1 —

and was monitored by AFTAC, Patrick AFB, FL 32925, under Contract No. F08606-77-C-000l. yax,z(x,) (02 + 72)2

{3—x,[x,+2(x,)]} (A2.52)

Acknowledgement The research work presented in this paper was mitiated while the author was visiting the Applied Seismology Group, Uncoln Laboratory, Massachusetts Institute of Technology. I would like to express my thanks to M. Chinnery and R.T. Lacoss for their valuable discussions and suggestions. I also wish to thank R. Needham and D. Newton for data and computer support. Finally I would like to thank E.S. Husebye and F. Ringdal at NORSAR for many valuable discussions in connection with this work. This research was supported in part by the Advanced Research Projects Agency of the U.S. Department of Defense,

References Aki, K., 1965. Maximum likelihood estimate of b in the formula log N a — bM and its confidence limits. Bull. Earthquake Res. Inst., Tokyo Univ., 43: 237—239. Christoffersson, A., 1970. The one component model with incomplete data. Department of Statistics, University of

=

Uppsala, 138 pp. Herrin, E. and Tucker, W., 1972. On the estimation of bodywave magnitude. Technical Report to AFOSR, Dallas

Geophysical Southern Methodist University, Dallas, TX, 44Laboratory, pp.

Kelly, E.J. and Lacoss, R.T., 1969. Estimation of seismicity and network detection capability. Tech. Note 1969-4 1,

Lincoln Laboratory, Massachusetts Institute of Technol. ogy, Cambridge, MA. Ringdal, F., 1975. the estimation of 65: seismic detection thresholds. Bull.On Seismol. Soc. Am., 1631—1642.

Ringdal, F., 1976. Maximum likelihood estimation of seismic magnitude. Bull. Seismol. Soc. Am., 66: 789—802.