Troyan V., Kiselev Y. Statistical Methods of Geophysical Data Processing

Подождите немного. Документ загружается.

Chapter 6

Methods for parameter estimation of

geophysical objects

At a choice of method for the parameters estimation, it is necessary to implement

study of properties of the estimated values. The main properties of the estimate were

considered in § 2.2, they are: the consistence, unbiasedness, eﬃciency (minimum

of variance), robustness and normality. Some of these properties can be executed

asymptotically, i.e. at a great number of measurements, for example, asymptotic

unbiasedness, asymptotic consistence, asymptotic normality. Let’s consider meth-

ods of a point estimation most frequently used at processing of real geophysical

data.

6.1 The Method of Moments

The method of moments is reduced to following. Let’s we have a repeated sample

, x

, . . . , x

of a random variable ξ with the likelihood function

L(x

, x

, . . . , x

; θ

, . . . , θ

) =

i=1

p(x

; θ

, . . . , θ

). (6.1)

To ﬁnd S sample moments

˜α

i=1

, s = 1, 2, . . .

If α

= M[x

] exists, then at n → ∞ the law of large numbers is valid, and we can

write

P [|˜α

− α

| > ε] → 0 if n → ∞,

i.e. at a large sample size ˜α

tends to α

with probability close to 1. Therefore it

is natural to construct the equations

∞

−∞

p(x; θ

, . . . , θ

)dx = α

, s = 1, 2, . . . , (6.2)

which allow to ﬁnd the estimates of the desired parameters θ

, . . . , θ

. At n → ∞

these estimates are asymptotically unbiased and asymptotically normal.

163

164 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

6.2 Maximum Likelihood Method

For any sample (x

, x

, . . . , x

) the likelihood function (6.1) is a function of an

unknown parameter vector θ. The maximum likelihood estimate of θ is a quantity

θθ, which corresponds to a maximum of the likelihood function ln L(x

, x

, . . . , x

; θ).

We can obtain the vector of estimates

θθ by the solution of the system of S likelihood

equations:

∂ ln L(x

, . . . , x

; θ)

∂θ

= 0, s = 1, . . . , S, (6.3)

if the likelihood function is a diﬀerentiable function.

In spite of the fact that the maximum likelihood method leads to the more

complicate computation in comparison of the method of moments, the maximum

likelihood estimator is preferable (especially in a case of a small sample size) due

to the following theorems, which we shall give without the proof.

(1) If the eﬀective estimate exists, it will be the unique solution

θ of the system of

likelihood equations (6.3).

(2) If the suﬃcient estimate exists, each solution of the system of likelihood equa-

tions (6.3) will be a function of this estimate. At suﬃciently general conditions

of the regularity of the function ln L, that are usually executed in practice, the

maximum likelihood estimator has following properties.

(1) At a great enough number of observations the solution of the likelihood

equation

θ or the system of likelihood equations

θθ (in case of a parameter

vector) converges in probability to a true value θ or

θ, i.e. the maximum

likelihood method gives the consistent estimate.

(2) The estimate of the maximum likelihood method is an asymptotically ef-

ﬁcient estimate, i.e. at a great number of observations the Rao-Cramer

inequality transfers to the equality, thus the variance of the estimated pa-

rameter reaches the lower edge:

n→∞

= [I

(F )

(θ

)]

−1

in the case of one parameter θ,

θ)

n→∞

= [I

(F )

(θ

)]

−1

in the case of the vector parameter θ.

(3) The maximum likelihood estimate

θ is an asymptotic normal estimate, i.e.

√

−θ

θ) ∼ N(0, [I

(F )

(θ

)]

−1

)

if n → ∞.

(4) The maximum likelihood estimate is an asymptotic unbiased estimate, i.e.

] = θ

if n → ∞.

Methods for parameter estimation 165

For a ﬁnite sample size n there is only one case, when the maximum likelihood

estimator is optimal — it is the exponential shape of the distribution function

of the general population. The maximum likelihood method is widely used for a

ﬁnding of the estimates of parameters of geophysical objects in the problems of the

quantitative interpretation. Thus, as a rule, the additive model of the experimental

data is used

u = {u

, u

, . . . , u

}, uu = ff (

θ) +

ε, (6.4)

where f(θ) is the model of the geophysical ﬁeld under investigation, which depends

on the vector of desired parameters, ε is a random error of measurement. The

normality hypothesis of the random component ε (εε ∈ N (0, R

)) allows, in the case

of an additive model (6.4), to write a logarithm of the likelihood function, without

the terms which do not contain θθ, in the form

(

u,θθ) = −

(uu −ff(

θ))

−1

(uu −ff(

θ)). (6.5)

It is necessary to note, that in case of the normal distribution, the estimates of the

maximum likelihood method are equivalent to the estimates of the weighted least

squares method where as the weights appear the elements of an inverse covariance

matrix R

−1

For the noncorrelated random component εε and observations with an equal

accuracy, the logarithm of the likelihood function (6.5) we can rewrite in the form

(uu, θ) = −

2σ

i=1

− f

(θ))

. (6.6)

The maximum likelihood criterion in this case turn to the classical least square

method. Substituting the relation (6.6) to the system of equations (6.3), we obtain

(

u −

f(θ))

−1

∂f

∂θ

= 0, s = 1, 2, . . . , S. (6.7)

6.3 The Newton–Le Cam Method

Let’s consider an approximate method of the solution of the system of equations

(6.7), which one is update of the Newton iterative method of the solution of the

system of algebraic equations. The idea of the method is reduced to following. The

initial parameter vector θ

(0)

is picked from a priori data. Further, the logarithm of

the likelihood function in the vicinity of the point

(0)

is expanded in the Taylor

series with three ﬁrst terms

(θθ) ≈ l

(θ

(0)

) + ∆θθ

d −

∆θ

C∆θ, (6.8)

where ∆θθ

= (θ −θ

)

is a diﬀerence of two vectors (row vector), dd is a derivative

on parameters (column vector), C is a matrix of the mathematical expectations of

the second derivatives

∂l

(θ

θ)

∂θ



θ=θ

(0)

, c

= −

∂

(θ

θ)

∂θ



θ=θ

(0)

166 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Substituting the expression (6.8) to the system of equations (6.3), after simple

transformations we obtain the system of linear likelihood equations in a form

C∆θθ = d. (6.9)

This system of the linear equations can be solved, for example, by the Gauss method,

if C is an improper matrix. In the case of a singularity of this matrix, it is necessary

to use the singular analysis or regularization methods(see § 6.12).

By solution of the system (6.9) we obtain the estimate of the vector ∆θθ,

∆

(1)

= C

−1

and a corrected value of the parameter vector

(1)

(0)

+ ∆

θθ

(1)

The value of the superscript is a number of the iteration. Further, the value of the

vector

(1)

becomes as an initial vector. We calculate the tare vector ∆

(2)

and

the corrected vector

(2)

The iterative process is prolonged until the relative variation of the vector ∆

θθ

(n)

becomes less than a priori threshold value β:



∆

(n)



≤ β, s = 1, . . . , S,

usually β is order of 10

−2

÷10

−3

. The rate of convergence depends on the function

f(θ

θ) and on the initial vector

(0)

In the considered algorithm the elements of the matrix C from (6.8) contain the

second derivative of the function

θ) on components of vector

θ. In a practical

sense it is more preferable instead of the matrix C to use the matrix

C = M[C],

which obtains by averaging of the matrix C over the elements of the random vector

εε with the next assumption M[ε] = 0. In this case the elements of matrix

C do not

contain the second derivatives of the function f(θθ). Let’s write the expression for

the elements of matrix

C belonging to n-th iteration

˜c

(n)

= M



−(

u −

f(θ))

−1

∂

∂θ



θ=

θθ

(n)

∂

∂θ



θ=

(n)

−1

∂f

∂θ



θθ=

(n)



= −M





−1

∂

∂θ

∂

∂θ



θ=

(n)

−1

∂f

∂θ



θ=

(n)

∂

∂θ

−1

∂

∂θ



θ=

(n)

Methods for parameter estimation 167

Let us note, that at the mathematical derivation of the components of matrix

C we

take into account the next equality M[ε] = 0 for the random vector ε.

The matrix

C does not depend on experimental data and completely deﬁned by

the model of the determinate and random components.

Using the model of experimental data and the formula (6.8), to ﬁnd an explicit

form of vector

= (uu −ff(

θ))

−1

∂f

∂θ



θ=

(n)

If the function

f(θ) is a linear one relatively to the desired parameters and we

work under the model

u = ψθ +

ε, (6.10)

where ψ is a plan matrix (or structure matrix) with dimension [n × S]:

ψ =



. . . ψ

. . . . . . . . . . . .

. . . ψ



then the system of equations (6.7) is a linear one relative to θ and it can be written

in a matrix shape

(ψ

−1

ψ)

θ = ψ

−1

It is easy to ﬁnd the maximum likelihood estimate, if the matrix ψ

−1

ψ is non-

singular:

θ = (ψ

−1

ψ)

−1

uu. (6.11)

In the case of singularity it needs to use the singular analysis or regularization (see

Sec. 6.12).

We have demonstrated, that in the case of the normal distribution and at pres-

ence of correlations, the maximum likelihood method is equivalent to the weighted

least squares method, and for independent components of the random vector ε

ε it is

equivalent to the least square method.

6.4 The Least Squares Method

The maximum likelihood method (MLM), as it was mentioned in § 6.2, at a normal

distribution of the random component of the model is equivalent to the least squares

method (LSM). The computing circuit of LSM is one of the simplest, therefore the

method is used for ﬁnding the estimates even in the case when the distribution of

experimental data diﬀers from the normal distribution. The estimates, obtained by

LSM in this case have smaller accuracy, but it is compensated by the simplicity of

their deriving.

168 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Let us consider a special case of LSM which has a great practical importance.

This is the case of the linear model relatively to the vector of desired parameters

(6.10). The desired estimates is found by minimization of the quadratic form

λ(

θ) = (

u − ψ

θ)

W (u

u − ψθ

θ),

θ = arg min λ(θ

θ),

where W is a weight matrix. The necessary condition consists in making to zero

the partial derivatives of the quadratic form on variables θ

, θ

, . . . , θ

∂λ

∂θ

= 0, s = 1, 2, . . . , S. (6.12)

The obtained system of equations is called the system of normal equations. Under

accepted assumptions, the system of equations (6.12) it is possible to write as

(ψ

W ψ)

θ = ψ

Wuu.

Using this relation it is easy to obtain LSM estimates, if the matrix ψ

W ψ is not

singular one:

θ = (ψ

W ψ)

−1

u. (6.13)

The estimate (6.13) coincides with MLM estimate (6.11) if the weight W matrix is

equal to the inverse covariance matrix R

−1

. In the case of an uncorrelated vector

of experimental data u, the main diagonal of the matrix R

−1

contains the inverse

variances σ

of the measurements u

. The analogy with MLM to clariﬁes a sense

of the weight matrix W , which appears formally in LSM. In a case of observations

with the equal accuracy

W = σ

−2

the estimate (6.13) can be written as

θθ = (ψ

ψ)

−1

ψu. (6.14)

The covariance matrix

θθ

= M

(

θ −θ

θ)(

θ −θ

θ)

= (ψ

ψ)

−1

ψ(ψ

ψ)

−1

is obtained from (6.14) by substituting R

= σ

I. Finally, we obtain

= σ

(ψ

ψ)

−1

The LSM estimates in the case of the linear model have the next optimal prop-

erties which do not depend on the distribution law of

ε.

(1) Unbiasness of

θ. Taking into account that M[u] = ψθθ, the mathematical ex-

pectation of the estimate

θ] = (ψ

ψ)

−1

M[u]

is equal to the value of the desired vector parameter θθ

θθ] = (ψ

ψ)

−1

ψθ = θ.

Methods for parameter estimation 169

(2) Eﬃciency of

θ. Among the class of unbiased estimates

∗

of the parameter

vector θθ, which are a linear combination of the data uu, LSM estimate

θ has a

minimum variance

(

θ) ≤ R

(

∗

)

for arbitrary s = 1, . . . , S.

(3) An estimate ˆσ

using residual sum of squares is given by the formula:

ˆσ

= (u

u − ψ

θ)

(

u − ψ

θ)/(n − S).

6.5 LSM — Nonlinear Model

The LSM estimate for the linear model has a simple computing circuit, but it also

can be used as a basis for the general case of the nonlinear model:

uu = ϕϕ(θ) + ε. (6.15)

Let’s represent ϕϕ(θ) as an expansion in Taylor series in vicinity of some initial

vector

and keeping two terms of the expansion, we rewrite the model (6.15) as:

˜u

u = ψ∆

θ + εε, (6.16)

where

˜u

u = u − ϕ(

), kψ

k =

∂ϕ

∂θ



θ=

, ∆θ

(0)

= θθ −θθ

The model (6.16) is equivalent to the model (6.10). The applying of the weighted

LSM, gives the estimate

∆

θθ = (ψ

W ψ)

−1

We repeat this procedure, but instead of θ

we use θθ

= θ

+ ∆

θ.

The convergence of the iterative process depends on a type of the function

ϕ(θθ)

and depends on a choice of weights. The iterative process stops, when the relative

variation of the arguments on the next step k becomes less than a given threshold

value β



∆θ

(k)



≤ β.

Usually the value of β is equal to 10

−2

– 10

−3

170 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

6.6 LSM — Orthogonal Polynomials (Chebyshev Polynomials)

The inversion of the matrix ψ

ψ for a derivation of the estimate (6.14) is simpliﬁed

essentially, if as the determinate function the system of the normalized orthogonal

polynomials is used. In this case i-th row of matrix ψ can be represented as

= kϕ

), ϕ

), . . . , ϕ

S−1

)k,

where ϕ

) are orthogonal on a set of x:

i=1

)ϕ

) =







i=1

), if s = s

0, if s 6= s

(6.17)

From the condition (6.17) it yields that the matrix ψ

ψ is a diagonal one

ψ =



i=1

) 0 . . . 0

i=1

) . . . 0

. . . . . . . . . . . .

0 0 . . .

i=1

S−1

)



and it is easy to be inverted.

We will consider the Chebyshev polynomials. The scheme of their construct-

ing looks like the next. Let us assume that ϕ

) = 1. The function ϕ

) is

determined as

) = x

+ b

). (6.18)

To ﬁnd the coeﬃcient b

using the orthogonality condition

i=1

)ϕ

) = 0, b

i=1

)

i=1

)

. (6.19)

Substituting b

from the formula (6.19) to the right hand side of the equality (6.18),

we obtain

) = x

−

i=1

)

i=1

)

Using the expression

) = x

+ b

) + b

Methods for parameter estimation 171

to ﬁnd ϕ

). Using the orthogonality condition to ﬁnd b

and b

i=1

)ϕ

) = 0, b

= −

i=1

)

i=1

)

i=1

)ϕ

) = 0, b

= −

i=1

)

i=1

)

Finally, for ϕ

) we obtain

) = x

−

i=1

)

i=1

)

) −

i=1

)

i=1

)

By analogy we can determine ϕ

) and so on. For an arbitrary s-th function we

have

) = x

−

i=1

s−1

)

i=1

s−1

)

s−1

) − ··· −

i=1

)

i=1

)

). (6.20)

It is easy to show the validity of a such expression. Let’s write the function

) in a following form:

) = x

+ b

s−1

) + ···+ b

After multiplying by ϕ

), s

< s, summing over x

, . . . , x

and using the orthog-

onality condition (6.17), we obtain

i=1

) + b

i=1

) = 0, b

= −

i=1

)

i=1

)

that demonstrates validity of the equality (6.20). For a case, when the abscissas

, . . . , x

are preselected with an uniform sample, the recursion formula

s+1

) =

−

i=1

) −

− i

)

4(4i

− 1)

s−1

)

for the calculation of the function ϕ

s+1

) on known ϕ

) and ϕ

s−1

) was

obtained by Chebyshev.

We have considered alternatives of a representation of the matrix ψ by orthog-

onal polynomials for the homogeneous process.

172 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

The estimate of the coeﬃcients of orthogonal polynomials by the LSM method,

using (6.13) and orthogonality conditions (6.17) under W = 1,

j=1

= 1 has a

form

j=1

, s = 0, 1, 2, . . . , S − 1,

θ = ψ

In this case, the covariance matrix of the estimated parameters is diagonal one

= σ

and the standard deviation gives by the formula

ˆσ

(

u − ψ

θ)

(u − ψ

θ)

n − S

6.7 LSM — Case of Linear Constraints

At the solution of applied problems a necessity frequently appears to ﬁnd the esti-

mates of desired parameters at presence of the linear constraints. It is one of the

alternatives to incorporate a priori information.

Let’s a vector of desired parameters θ satisﬁes the system of Q equations

Aθ = V , (6.21)

where Q < S. The matrix A has a dimension [Q × S] and the vector V

V has a length

Q. To determine the desired parameters the problem can be stated as a problem on

the conditional extremum, i.e. to minimize of λ(θ) under the condition (6.21). For

the solution of this problem, as usual, we will use the Lagrange multipliers method.

To build functional

Φ(θ

θ) = (u

u − ψθ

θ)

(

u − ψ

θ) − 2

(Aθ

θ −

V ),

where a component of the vector λ is the Lagrange multiplier. The length of λλ

is equal to the length of V

V , i.e. the lengths of the vector are equal to a number

of linear constraints. As the estimate of the parameters we will choose the values,

under which the functional Φ(θ

θ) reaches an own minimum. The necessary condition

of a minimum consists in the equalling to zero of the ﬁrst derivative on components

of the vector

θ:

∂Φ(

θ)

∂θ

= 0, s = 1, . . . , S. (6.22)

After the solution of the system of linear equations (6.22) relatively to θ

, we obtain

the next estimate

θ, expressed using Lagrange multipliers:

+ λ

A(ψ

ψ)

−1

, (6.23)