Troyan V., Kiselev Y. Statistical Methods of Geophysical Data Processing
Подождите немного. Документ загружается.
Statistical criteria for choice of model 203
of qualitative interpretation is applied the criterion of a posteriori probability ratio
for the model (7.1) at given parameter vectors
ˆ
θ
0
and
ˆ
θ
1
.
Let us consider the hypotheses H
0
and H
1
with given a priori probabilities p(0)
and p(1) respectively. Using the Bayes theorem for the calculation of a posteriori
probabilities to find the criterion function
λ = ln
P (1/uu)
P (0/uu)
= ln
P (1)p(u/1)
P (0)p(u/0)
.
An application of the criterion consists in a test of inequalities
if λ < 0, then H
0
does not contradict u,
if λ ≥ 0, then H
0
is rejected.
The criterion function λ is a random one, because it depends on a random measured
data. Hence, a most completely description of λ is reduced to the determination
of its distribution function. For a simplicity of the analysis, we assume that ε
0
∈
N(0, R),
ε
ε
1
∈ N(0, R), where R is the covariance matrix of the random component,
then λ will be written as
λ =
1
2
(u − f(θ
0
))
T
R
−1
(uu −
f
f(θ
0
)) −
1
2
(u −ff (θθ
1
))
T
R
−1
(u − f(θ
1
)). (7.4)
To find the conditional mathematical expectations of λ
0
, λ
1
and the conditional
variances of σ
2
λ
0
, σ
2
λ
1
, using the following models
u = f(θ
0
) +
ε
ε,
u
u =
f
f(θ
1
) + εε,
hλ
0
i =
1
2
hε
T
R
−1
εi −
1
2
(f(θθ
0
) − f (θ
1
))
T
R
−1
(ff(
θ
θ
0
) − f(θθ
1
))
−
1
2
hε
T
R
−1
εi − h(f(θ
0
) −ff(
θ
θ
1
))
T
R
−1
εi
= −
1
2
(f(θ
0
) − f(θ
1
))
T
R
−1
(ff(θθ
0
) − f (θ
1
)) = −hλ
1
i, (7.5)
σ
2
λ
0
= h(λ
0
− hλ
0
i)
2
i = h[(−
1
2
(f(θ
0
) −
f
f(θ
1
))
T
R
−1
× (ff(θθ
0
) − f (θ
1
))) − (f (θ
0
) − f(θ
1
))
T
R
−1
ε
+
1
2
(f(θ
0
) −
f
f(θ
1
))
T
R
−1
(f(θθ
0
) − f (θ
1
))]
2
i
= (ff(θθ
0
) − f (θ
1
))
T
R
−1
h
ε
εε
T
iR
−1
(ff(
θ
θ
0
) − f (θθ
1
))
= (ff(θθ
0
) − f (θ
1
))
T
R
−1
(
f
f(
θ
θ
0
) − f(θθ
1
)) = σ
2
λ
1
. (7.6)
From an analysis of the formula (7.4) leads that λ is a linear function of the measured
data u. Hence, at mentioned above assumptions, λ has the normal distribution
p(λ/0) =
1
√
2πσ
λ
0
exp
−
(λ − hλ
0
i)
2
2σ
2
λ
0
,
for the null hypothesis H
0
and
p(λ/1) =
1
√
2πσ
λ
1
exp
−
(λ − hλ
1
i)
2
2σ
2
λ
1
for the hypothesis H
1
. In accordance with the criterion, if λ belongs to the interval
[−∞, 0], then the hypothesis H
0
is not rejected. If λ belongs to the interval [0, ∞]
it is not rejected the alternative hypothesis H
1
. The next four special cases can
appear (Fig. 7.1):
204 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
(1) the hypothesis H
0
is true and it is not rejected;
(2) the hypothesis H
0
is true, but it is rejected (error of the first kind, alpha error);
(3) the hypothesis H
0
is false and it is rejected;
(4) the hypothesis H
0
is false, but it is not rejected (error of the second kind, beta
error).
Fig. 7.1 The density functions p(λ/0) (hypothesis H
0
) and p(λ/1) (hypothesis H
1
) and decision
errors.
Let denote a probability of the error of the first kind as P
01
and a probability
of the error of the second kind as P
10
. To calculate these probabilities:
P
01
=
∞
Z
0
p(λ/0)dλ =
1
√
2πσ
λ
0
∞
Z
0
exp
−
(λ − hλ
0
i)
2
2σ
2
λ
0
dλ
=
1 − Φ
−
hλ
0
i
σ
λ
0
,
P
10
=
0
Z
−∞
p(λ/1)dλ =
1
√
2πσ
λ
1
0
Z
−∞
exp
−
(λ − hλ
1
i)
2
2σ
2
λ
1
dλ
= Φ
−
hλ
1
i
σ
λ
1
,
Φ(x) =
1
√
2π
x
Z
−∞
exp{−t
2
/2}dt, Φ(−x) = 1 − Φ(x),
where Φ(x) is the normal density function. To estimate a decision quality, proba-
bility of the total error is introduced
P
e
= P (0)P
01
+ P (1)P
10
.
For the considered special case of the equality of a priori probabilities: P (0) =
= P (1) = 1/2 by hλ
0
i = −hλ
1
i, σ
λ
0
= σ
λ
1
we obtain
P
e
=
1
2
1 − Φ
−
hλ
0
i
σ
λ
0
+
1
2
Φ
hλ
0
i
σ
λ
0
= Φ
hλ
0
i
σ
λ
0
.
Substituting (7.5) (λ
0
) and (7.6) (σ
λ
0
) in above expression, we obtain
P
e
= Φ(−α/2) = 1 − Φ(α/2) (7.7)
Statistical criteria for choice of model 205
α = [(
f
f(
θ
θ
0
) −f
f
f(
θ
θ
1
))
T
R
−1
(
f
f(
θ
θ
0
) −
f
f(
θ
θ
1
))]
1/2
. (7.8)
At the calculation we should set an allowed value of the error probability P
e
(unreliability of the problem) and to find a value of α
0
(Fig. 7.2). If we know this
Fig. 7.2 Graphic representation of the function Φ(−α/2).
value, then in all cases of the validity of the inequality
[(
f
f(
θ
θ
0
) −f
f
f(θ
θ
θ
1
))
T
R
−1
(
f
f(
θ
θ
0
) −f
f
f(
θ
θ
1
))]
1/2
≥ α
0
,
the unreliability can not be greater than P
e
. So, a value of α
0
determines a threshold
condition for the decision with the given efficiency.
Specifically, at P
e
= 0.05, α
0
= 3.3 the threshold condition looks like
α = [(
f
f(
θ
θ
0
) −f
f
f(
θ
θ
1
))
T
R
−1
(
f
f(θ
θ
θ
0
) −
f
f(
θ
θ
1
))]
1/2
≥ α
0
= 3.3.
The inequalities obtained above allow us to determine the analytical threshold con-
ditions between the parameters of the considered problem in the case if the unreli-
ability is smaller than 0.05.
Example 7.1. Let us consider an unreliability of the extraction of the single regular
wave on a noise background (Fig. 7.3):
u
ki
=
either Aϕ(t
i
− τ − k∆xγ) + ε
ki
,
either ε
ki
,
where k = 1, 2, . . . , K (K is a number of geophones); ∆x is a distance between
geophones; A, τ, γ are the amplitude, the arrival time, the apparent velocity of
a seismic wave respectively. The unreliability P
e
is determined by the expression
Fig. 7.3 The regular wave.
(7.7), in addition, α has the next form:
α
2
=
A
2
σ
2
X
k
X
i
X
i
0
r
−1
ii
0
ϕ(t
i
− τ − k∆xγ)ϕ(t
i
0
− τ − k∆xγ),
206 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
where r
−1
ii
0
are the elements of the inverse covariance noise matrix, σ
2
is the noise
variance. In the specific case of the uncorrelated noise component we obtain
α
2
=
A
2
σ
2
n
X
i=1
K
X
k=1
ϕ
2
(t
i
− τ − k∆xγ).
7.3 The Signal Resolution Problem
Let us consider the model for the signal resolution
u
ki
=
either A
0
ϕ(t
i
− τ
0
− k∆xγ
0
) + ε
ki
,
or A
1
ϕ(t
i
− τ
1
− k∆xγ
1
) + A
2
ϕ(t
i
− τ
2
− k∆xγ
2
) + ε
ki
.
The unreliability of the waves separation P
e
(or probability of erroneous decision)
is described by the formula (7.7), where
α
2
=
1
σ
2
X
k
X
i
X
i
0
r
−1
ii
0
(A
0
ϕ(t
i
− τ
0
− k∆xγ
0
) − A
1
ϕ(t
i
− τ
1
− k∆xγ
1
)
− A
2
ϕ(t
i
− τ
2
− k∆xγ
2
))(A
0
ϕ(t
i
0
− τ
0
− k∆xγ
0
)
− A
1
ϕ(t
i
0
− τ
1
− k∆xγ
1
) − A
2
ϕ(t
i
0
− τ
2
− k∆xγ
2
)).
For the case of the uncorrelated random component, at the symmetric arrangement
of the waves f
1
and f
2
with respect to f
0
, i.e.
τ
1
= τ
0
− ∆τ, γ
1
= γ
0
− ∆γ,
τ
2
= τ
0
+ ∆τ, γ
2
= γ
0
+ ∆γ,
and under the condition A
0
= A
1
/2 = A
2
/2 we have
α
2
=
A
2
0
σ
2
X
k
X
i
(ϕ(t
i
− τ
0
− k∆xγ
0
) −
1
2
ϕ(t
i
− (τ
0
− ∆τ)
− k∆x(γ
0
− ∆γ)) −
1
2
ϕ(t
i
− (τ
0
+ ∆τ) −k∆x(γ
0
+ ∆γ)))
2
.
The analysis of α
2
can lead that the probability of erroneous decision (unreliability)
depends on: a choice of the observation scheme (∆x and K), the signal-to-noise ratio
(A/σ), the signal shape ϕ(t), the kinematic parameters of waves (τ
0
, γ
0
, τ
1
, γ
1
, τ
2
,
γ
2
). Studying a dependence of α
2
on these parameters it is possible to solve the
next problems.
(1) To determine source-to-receiver offset distance, provided the given kinematic
resolution of waves with the assigned probability error P
e
(at the fixed signal-
to-noise ratio).
(2) To estimate the minimum acceptable signal-to-noise ratio, which gives a needed
resolution of the waves using the kinematic parameters with the assigned unre-
liability for the chosen observation scheme.
Statistical criteria for choice of model 207
Example 7.2. We consider the function α
2
(∆τ) corresponds to the signals ϕ(t) =
exp{−βt}cos(Ωt) under the condition γ
1
= γ
2
= γ
0
, i.e. resolution is studied only
for the time differences τ
1
= τ
0
−∆τ, τ
2
= τ
0
+ ∆τ. An example of a such function
is represented on Fig. 7.4. In the initial part of the function there is a maximum
at the time τ ∼ T /2 and after several oscillations the curve tends to an asymptote,
which corresponds to a waves reception without an interference.
Fig. 7.4 Graphic representation of the function α
2
.
Example 7.3. Let us consider an example of the resolution of magnetized bodies.
We will find the threshold condition for the discrimination of the vertical thin layer,
with the infinite embedded lower edge and located under the coordinate origin
transversely to the observation axis x, and two similar layers located symmetrically
relative to the coordinate origin at the distances ±ξ. The magnetization I of the
single object is equal to the doubled magnetization of each of the paired layer. So a
total magnetization of the pair of layers is equal to the magnetization of the single
layer. If the random component is uncorrelated at the observation points, then
f
0k
=
2I
h
1
1 + (k∆x/h)
2
, f
1k
=
I
h
1
1 + (k∆x/h − ξ/h)
2
,
f
2k
=
I
h
1
1 + (k∆x/h + ξ/h)
2
. (7.9)
An example of the magnetic fields produced by the considered objects are repre-
sented on Fig. 7.5. The threshold condition can be written as
K
X
k=−K
2
1 +
k∆x
h
−
ξ
h
2
−1
+
2
1 +
k∆x
h
+
ξ
h
2
−1
−
1 +
k∆x
h
2
−1
2
≥
σh
2I
2
α
2
0
.
The obtained inequality determines the minimum value of ξ, which still makes
possible to distinguish the interference field of the pairs of layers from the field of
the single layer.
208 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
Fig. 7.5 The magnetic field produced by the vertical layers (I = 1, ξ = 2.5, h = 1.0): (a), (c) —
without the noise; (b), (d) — with uncorrelated Gaussian noise N (0, 0.25); (a), (c) — one layer;
(b), (d) — two layers.
7.4 Information Criterion for the Choice of the Model
Let us consider the information criterion (Akaike, 1974). The application of the
maximum likelihood method for the problem of the choice of the model at the given
set of allowable models with a different dimension of the parameter vectors, leads to
the choice of the model with greater dimension of the parameter vector. Therefore
this method can not be used for the intuitive chose of the “right” model. However,
the analysis of the maximum likelihood estimates results in that, taking into account
the asymptotic properties of the maximum likelihood estimate, the likelihood func-
tion is possible to consider as the keen value relatively to the estimated parameters
in a vicinity of their true quantities. This property can be used for the development
of the fitting criterion of the model with the probability structure described by the
density function p
1
(u,
θ
θ) with respect to the probability structure described by the
density function p
0
(u). In the case of the classical maximum likelihood method,
the estimation is implemented in the framework of one family of distributions. The
information criterion allows belonging of p
1
(
u
u,θθ) and p
2
(u) to the family with vari-
ous dimensions of the desired parameter vectors or to various families. At this case
the problem of the choice of the model is stated as the estimation problem and we
are free from the subjective decisions.
Let us consider the construction of the information criterion. Let the values of
the observed geophysical field u
1
, u
2
, . . . , u
n
are independent random values with
a priori unknown density function p
0
(
u
u). The parametric family p
1
(u, θ) has an
unknown parameter vector, moreover, we do not know both the components of the
vector θ, and its dimension. The estimate of the mean value of the logarithm of
Statistical criteria for choice of model 209
likelihood function has the form
1
n
n
X
i=1
ln p
1
(u
i
, θ). (7.10)
Passing on to the limit in the formula (7.10) under n → ∞, we obtain the mean
value of the logarithmic likelihood function
B[p
0
(u), p
1
(
u
u,θθ)] = lim
1
n
n
X
i=1
ln p
1
(u
i
,
θ
θ)
=
∞
Z
−∞
p
0
(
u
u) ln p
1
(u
u
u,
θ
θ)du
1
. . . du
n
.
Taking into account the asymptotic efficiency of the maximum likelihood esti-
mate, we can conclude that the mean value of the logarithmic likelihood function
B[p
0
(u
u
u), p
1
(
u
u,
θ
θ)] is a most sensitive measure with respect to the small deviations
of p
1
(u
u
u,θ
θ
θ) from p
0
(u
u
u). The difference
I
k
(p
0
(
u
u), p
1
(uu, θ)) = B[p
0
(u), p
0
(
u
u)] − B[p
0
(u), p
1
(
u
u,θθ)]
is known as the Kullback information (Kullback, 1978), which is a measure of the
deviation of the actual density function p
0
(u) from the expected one p
1
(u,
θ
θ). If
p
1
(u, θ) everywhere is equal to p
0
(u), then the information I
k
is equal to zero.
Let us consider a special case. Let p
0
(u) = p
0
(
u
u,θθ
0
), where θ
0
is the actual
parameter vector. The mean value of the logarithmic likelihood function and the
Kullback information depend only on the parameter vectors θθ
0
and θθ:
I
k
(θ
0
,θθ) = B(θθ
0
,θθ
0
) − B(θ
0
, θ),
where
B(θθ
0
,θθ
0
) =
∞
Z
−∞
p
0
(u),
θ
θ) ln p(uu), θθ
0
)du
1
. . . du
n
.
If the vector θ
0
is sufficiently close to the vector θθ (θθ = θ
0
+ ∆θθ), then it is possible
the next approximation
I
k
(θ
0
, θ
0
+ ∆θ) =
1
2
∆θ
T
I
F
∆θ, (7.11)
where
kI
F
k
ss
0
=
∂ ln p(u, θ)
∂θ
s
∂ ln p(uu,θθ)
∂θ
s
the elements of the Fisher matrix. After the analysis of the formula (7.11) we can
conclude, that if the estimate of the parameter
ˆ
θ is located in a small neighborhood
of the parameter θ
0
, then the deviation of the density p
1
(uu, θ) from the actual
density p
0
(u,
θ
θ
0
) is measured by the mean square deviation
1
2
(
ˆ
θ
θ
θ −
θ
θ
0
)
T
I
F
(
ˆ
θ
θ
θ −θ
θ
θ
0
).
210 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
If θ
θ
θ is closed to θ
θ
θ
0
in the maximum point of B(θ
θ
θ
0
,θ
θ
θ):
ˆ
θ = arg max B(θ
0
, θ),
it is possible to proof, that under some regularity condition the maximum likelihood
estimate
ˆ
θθ has an asymptotically normal distribution, i.e.
√
n(
ˆ
θ −
θ
θ) ∈ N(0, I
−1
F
), (7.12)
n(
ˆ
θ − θ)
T
I
F
(
ˆ
θ − θ) ∈ χ
2
Sm
, (7.13)
where χ
2
Sm
is χ
2
-distribution with a number of degrees of freedom equated to a
length the vectors
ˆ
θ and θ, which belongs to the same parametric space S
m
. But
the vector θ
0
has the length S
m
0
, which is not equal to the length of the vectors
ˆ
θθ
and θ, in general.
Let us evaluate the mean Kullback information for the maximum likelihood
estimate
ˆ
θθ. To rewrite the formula (7.11) in the form
I
k
(θθ
0
,
ˆ
θθ) =
1
2
(
ˆ
θ − θ)
T
I
F
(
ˆ
θ − θ). (7.14)
To use the condition (7.13), we multiply the left hand side and the right hand
side of (7.14) by 2n. After simple reduction we will have
2nI
k
(θ
0
,
ˆ
θ) = n(
ˆ
θ −
θ
θ)
T
I
F
(
ˆ
θ
θ −
θ
θ) + n(
θ
θ −
θ
θ
0
)
T
I
F
(θ − θ
0
)
+ 2n(
ˆ
θ
θ −θ
θ
θ)
T
I
F
(
θ
θ −θ
θ
θ
0
). (7.15)
In accordance with the property (7.13), the first term from the right hand side of
the equality (7.15) belongs to the χ
2
S
m
distribution asymptotically, and the third
term has an asymptotically normal distribution. The calculation of the mean value
of the left hand side and the right hand side of the equation (7.15) on the ensemble
of observed data, we obtain
h2nI
k
(
θ
θ
0
,
ˆ
θ)i = n(θ − θ
0
)
T
I
F
(θθ − θ
0
) + S
m
.
To find n(
θ
θ −θ
0
)
T
I
F
(θθ −θθ
0
) we use the estimate (
θ
θ
0
,
θ
θ
0
) and (
θ
θ, θ
0
) in the form
of the sample mean of the logarithmic likelihood function. Substituting instead of
θθ its estimate
ˆ
θθ and to correct an arising bias by above way, we obtain:
n(
θ
θ −
θ
θ
0
)
T
I
F
(θ
θ
θ −θ
θ
θ
0
) = 2
n
X
i=1
ln p(u
i
,θ
θ
θ
0
) −
n
X
i=1
ln p(u
i
,
θ
θ) + S
m
. (7.16)
Substituting the expression (7.16) to the right hand side of the equality (7.15)
and taking into account that
I
k
(
θ
θ, θ
θ
θ) = B(θ
θ
θ
0
,θ
θ
θ
0
) − B(
θ
θ
0
,
θ
θ),
B(
θ
θ
0
,
θ
θ
0
) =
1
n
n
X
i=1
ln p(u
i
,
θ
θ
0
),
Statistical criteria for choice of model 211
we obtain the final expression for the function of the information criterion
IC(S
m
,
ˆ
θ) = −2hnB(
θ
θ
0
,
ˆ
θ
θ)i = −2 ln p(n,
ˆ
θ) + 2S
m
. (7.17)
The first term from the right hand side of this expression is equal to the max-
imum value of the likelihood function with an opposite sign. The second term is
equal to a doubled length of the vector θ. The function IC(S
m
,
ˆ
θ) can be considered
as discord measure of the identifiable model.
The criterion application is reduced to following. To choose the competing mod-
els m = 1, 2, . . . , M , with various lengths of the parameter vectors S
1
, S
2
, . . . , S
M
.
We calculate the estimates by the maximum likelihood method
ˆ
θ
m
and using the
formula (7.17) to find IC(S
m
,
ˆ
θ
θ). In accordance with the criterion, we choose the
model connected with a maximum value of IC, i.e.
S
m
= arg max
m
[IC(S
m
,
ˆ
θ
m
)]. (7.18)
The relevant advantage of the considered criterion is the actuation of the choice
of dimensionality of the model during the estimation. A quality of the choice is
estimated by an absolute value of the function IC in the maximum point on m.
The function IC is equal, within a constant factor, to the negative entropy. The
information criterion can be considered as a modification of the maximum entropy
method, which has a wide applying at the processing of geophysical data.
7.5 The Method of the Separation of Interfering Signals
The important step of the processing of geophysical data is the problem of the
separation of a geophysical field on its components. For example, one of the ob-
jectives of the analysis of a seismic field consists in the extraction of the reflected
waves which carry an information about parameters of the layered medium. In the
magnetometry the actual problem is the extraction of the fields produced by some
anomalies.
Let us consider the principles of the development of the iterative algorithm for
the separation of the interfering signals. We will find not only the signal parameters,
but to find also a number of the signals for the case of the representation of the
field model in the form
u
u
k
=
M
X
µ=1
f
f
f
µk
(
θ
θ
µ
) + ε
ε
ε
k
,
where θ
µ
is the unknown parameter vector for µ-th signal, εε
k
∈ N(0, R). For
example, in the case of seismic signals, we can write
f
µk
= kA
µ
ϕ(t
i
− τ
µ
− k∆xγ
µ
)k
n
i=1
.
To find the parameters we use the maximum likelihood method (see Sec. 6.2).
By analogy with the formula (6.5) we write the logarithmic likelihood function,
212 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
excluding the terms independent of θ,
l(u
k
, θ) =
1
2
X
k
(uu
k
−
M
X
µ=1
f
µk
(θθ
µ
))
T
R
−1
(u
k
−
M
X
µ=1
ff
µk
(θ
µ
)),
where θθ = kθ
µ
k
M
µ=1
is a full set of the vectors
θ
θ
µ
. Because u
k
does not depend on
the desired parameters θ
µ
, then it is enough to maximize the function
g(u
k
, θ) =
M
X
µ=1
X
k
uu
T
k
R
−1
f
f
µk
(θ
µ
) −
1
2
X
µ
0
6=µ
f
f
T
µ
0
k
(θθ)R
−1
f
µk
(
θ
θ)
−
1
2
f
T
µk
(
θ
θ
µ
)R
−1
f
µk
(
θ
θ
µ
).
, (7.19)
This function is a sufficient statistic and we will call it the sufficient reception
function. The first terms in the square brackets of the formula (7.19) is a correlation
between the signal f
µk
(θθ
µ
) of a number µ with measured data, the second term
determines a mutual correlation of the various signals and the third term is an energy
of µ-th signal. If the signals do not interact (by the reason of the orthogonality or
there is no overlapping between them in the space-time domain), then the second
term in square brackets is equal to zero and we can rewrite the formula (7.19) as
g(uu
k
,θθ) =
M
X
µ=1
X
k
uu
T
k
R
−1
f
f
µk
(θ
µ
) −
1
2
ff
T
µk
(θ
µ
)R
−1
ff
µk
(θ
µ
)
. (7.20)
As far as an each term of the sum over µ in the formula (7.20) depends only on
the desired parameter
θ
θ
µ
, the maximum of g(u
u
u
k
,θ
θ
θ) coincides with the maxima
of the terms of the sum over µ. It allows to find the parameters θ
θ
θ
µ
separately
(independently for each µ). In the case of correlated signals we should depart
from the above procedure and to vary the parameters simultaneously. In this case
the various iterative multi-step procedures can be offered with the variation only
parameters of one signal at each step.
For the algorithm development we rewrite the sufficient reception function (7.19)
and write down separately the terms corresponding to the m-th signal:
g(uu
k
,θθ)=
X
k
u
T
k
R
−1
f
mk
(θθ
m
) −
M
X
µ=1
µ6=m
X
k
f
µk
(
θ
θ
µ
)R
−1
f
mk
(θθ
m
)
−
1
2
X
k
f
T
mk
(θ
m
)R
−1
ff
mk
(
θ
θ
m
)
+
M
X
µ=1
µ6=m
X
k
u
T
k
R
−1
ff
µk
(
θ
θ
µ
)
−
1
2
M
X
µ=1
µ6=m
M
X
µ
0
=1
µ
0
6=m
X
k
f
f
T
µ
0
k
(
θ
θ
µ
0
)R
−1
f
f
f
T
µk
(θ
θ
θ
µ
)
−
1
2
M
X
µ=1
µ6=m
X
k
ff
T
µk
(θ
µ
)R
−1
ff
µk
(θ
µ
)
. (7.21)