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

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

Methods of transforms and data analysis 393

In practice the autocorrelation function R

(i) has a ﬁnite length. One can

suppose that R

(i) = 0, |k| > n outside the measuring span. But such values of

the autocorrelation function can be chosen so, that H is stable with respect to the

value R

(i) at |k| > n (Burg, 1972). Using the formulas (12.46), (12.48), we obtain

the system of equations:

∂H

∂R

(k)

= 0 =

−ω

∞

(z)

dω, |k| > n. (12.49)

We will written down the ﬁnal result without mathematical manipulations:

∞

(z) = E[h(z)h(z

−1

)]

−1

, (12.50)

where E and h(z) are the sum of predictive error squares (12.43) and z-transform

of the predictive error ﬁlter (12.44). The solution of the system of equations (12.44)

is considered in (Anderson, 1974; Smylie et al., 1973).

Consider a recurrent procedure for the solution of the system of equations (12.44)

to ﬁnd the predictive error ﬁlter 1, −h

, . . . , −h

(n = 0, 1, 2, . . . , N) in

the case of the time series x

, . . . , x

(m = 1, 2, . . . , M ) of a ﬁnite length.

We deﬁne the autocorrelation function taking into account a number of the terms

in considered time series:

M − m

M−m

k=1

k+m

(m = 0, 1, 2, . . . , M −1).

By analogy the sum of the squared diﬀerences of the predictive error ﬁlter E is

normalized.

Let’s note that the following identities

· 1 = E

, where E

k=1

in initial step n = 0 are satisﬁed. At the ﬁrst step (n = 1) the coeﬃcient h

of the

predictive error ﬁlter 1, −h

(where the ﬁrst number of the subscript is the step

number) is found by minimizing the sum of the prediction error

2(M − 1)

M−1

j=1

+ ˜e

which calculated by averaging of the predictions in the direct

= x

j+1

− h

and inverse

˜e

= x

− h

j+1

directions,

2(M − 1)

M−1

j=1



j+1

− h

)

+ (x

− h

j+1

)



. (12.51)

394 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

Considering E

as a function of h

, we ﬁnd h

as the solution of the equation

/dh

= 0

2(N − 1)

M−1

j=1

2[(x

j+1

− d

)(−x

) + (x

− d

j+1

)(−x

j+1

)]

2(N − 1)

M−1

j=1

[−(x

j+1

− x

j+1

) + h

j+1

+ x

)].

Thus,

= 2

M−1

j=1

j+1

M−1

j=1

+ x

j+1

Using expression (12.51), we ﬁnd E

. Let’s note, that R

it is possible to ﬁnd from

the ﬁrst equation of the system of equations



−h



Generally the average total error for the direct and inverse ﬁlters can be written

down as follows:

2(M − n)

M−n

j=1



−

k=1

j+k





j+n

−

k=1

j+nkk





The average error also can be written down as

2(M − n)

M−n

j=1

[(b

− h

)

+ (b

− h

)

where

= x

, b

= x

1+j

k=0

n−1,k

j+k

k=0

n−1,n−k

j+m−k

k=0

n−1,k

j+m−k

k=0

n−1,n−k

j+k

at that for the coeﬃcients b and b

and the coeﬃcients of the ﬁlter h the next

recurrent relations

= b

n−1,j

− h

n−1,n−1

n−1,j+1

= b

n−1,j+1

− h

n−1,n−1

n−1,j+1

= h

n−1,k

− h

n,m−k

, , k = 1, 2, . . . , n − 1,

= −1, h

= 0 (k > n)

Methods of transforms and data analysis 395

are valid. By solution of the equation for h

— dE

/dh

= 0, we obtain

= 2

M−n

j=1

M−n

j=1

+ b

). (12.52)

Let’s show that the value E

also satisﬁes to a recurrent relation. Let’s write

the system of equations determining the predictive error ﬁlter for the case of n = 3:



−h



. (12.53)

Substituting to the system of equation (12.53) the recurrent relation for h

, we

obtain















−h

n−1,1

−h

n−1,2



− h



−h

n−1,1

−h

n−1,2















n−1



− h



n−1



Generally the recurrent formula for E

can be written down as follows:

= E

n−1

(1 − h

Let’s note, that from equality (12.52) it follows |h

| < 1, therefore the inequalities

0 < E

< E

n−1

are valid.

This page intentionally left blankThis page intentionally left blank

Appendix A

Computer exercises

The computer exercises represented in the Appendix, make a basis of the computing

practical work which go with lectures “Statistical methods of processing of the

geophysical data” and “Geophysical inverse problems”, given to master students

and to post-graduate students in the Earth Physics Department of the Physical

Faculty of Saint-Petersburg State University within last ten years. The computer

exercises will allow the students to understand the essence and opportunities of the

considered methods of analysis and processing of the geophysical information more

deeply.

Basis of the computer exercises is the set of programs

, written in MATLAB

language. The programs contain a brief description of the input and output data,

demonstration example and recommendations for change input data.

The com-

puter exercises require the initial rudiments about MATLAB, such which allow to

change the input data in programs written in MATLAB language.

A.1 Statistical Methods

A.1.1 Examples of numerical simulation of random values

In Matlab package there are procedures for simulation of the random values. The

next script produces the samples belonging to uniform (at the interval (0, 1)) and

normal (N (0, 1)) distributions.

% Script r0a.m

% Examples of random values

% indd=1 : uniform distribution at interval (0,1)

% indd=2 : Gaussian distribution (N(0,1))

clear

% INPUT

The programs are placed in http://geo.phys.spbu.ru/

∼

vtroyan.

The programs together with demonstration examples was prepared in Geophysical Institute

of Lausanne University. Authors consider as the pleasant debt to express gratitude to professor

R.Olivier for creation of excellent conditions for work and a friendly climate.

397

398 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

indd=1;

m=25; % sample size

% END INPUT

n=1:m; % n=1,2,...,m

if indd==1

s=rand(1,m); % s : row vector

else

s=randn(length(n),1); % s : column vector; length(n)=m

s=s’ % s : row vector

end

rect=[0.1 0.1 .8 0.5]; % size of ﬁgure

axes(’position’,rect)

stem(n,s); % graphic representation

A.1.1.1 Exercises

(1) With the use of script r0a.m to implement a graphic representation of the

samples belonging to the uniform and normal distributions.

(2) By using the command line to implement graphic presentation of the samplings

belonging to the uniform and normal distributions with the use of plot com-

mand.

A.1.2 Construction of histogram

The histogram characterizes the probability properties of a random variable. The

next script (p0hist.m) constructs the histogram for a sampling belonging to the

normal distribution.

%Script p0hist.m

% Construction of histogram

% nr : sample size

clear

% INPUT

nr=500; % sample size

del=0.5; % input of bin for xcoordinate

% END INPUT

x=-5:del:5; % construction of the bins for x-coordinate

d=randn(1,nr); % d : a batch of numbers

%(row vector)

% hist(d,x);

[n,x]=hist(d,x);

bar(x,n./(nr.*del)); % graphic representation

% of the normalized histogram

Computer exercises 399

hold on

d2=d gauss(x,0,1); % normal probability density

% (row vector)

plot(x,d2,’g’) % graphic representation of % the normal density function

A.1.2.1 Exercises

(1) To compare the histograms obtained using the diﬀerent sample sizes (for exam-

ple, 100 – 500).

(2) To compare the histograms obtained for a diﬀerent number of the bins on the

x axis (for example, 4 – 20). Do not change the interval of the histogram

determination on the x axis.

A.1.3 Description of a random variable

A few probability distributions for the discrete random variables and the density

functions for the continuous random variables are represented below.

A.1.3.1 Beta distribution

Density function for the beta distribution reads as

(x|α, β) =

(

Γ(α+β)

Γ(α)Γ(β)

α−1

(1 − x)

β−1

, if 0 < x < 1,

0, if x ≤ 0 or x ≥ 1,

(A.1)

where α > 0, β > 0;

Mξ =

(α + β)

, Dξ =

αβ

(α + β)

(α + β + 1)

Graphic presentation of (A.1) is realized by the script d0beta.m. The function

d=d beta(x,al,bt) returns a row vector d with values of the density function for

the beta distribution.

A.1.3.2 Binomial distribution

Binomial distribution function is determined by the formula (1.38). Graphic pre-

sentation of the distribution function (1.38) is realized in the script d0bin.m. The

function d=d bin(n,p) returns row vector d — the probabilities of the binomial

distribution.

The script d0bin.m is represented below.

% Script d0bin.m

% Binomial distribution

% call function d=d bin(n,p)

% n : n=1,2,...; - the maximum value of x

% x=0,1,...,n

400 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

% p : 0 ¡ p ¡ 1

% d : row vector, contains the distribution function

clear

% INPUT

n=50;

p=0.7;

% END INPUT

d=d bin(n,p);

x=0:n;

% graphic representation

rect=[0.1 0.1 .8 0.5]; % size of the ﬁgure

axes(’position’,rect)

stem(x,d);

comp prob=sum(d) % computation of the total probability

math exp= n .* p % mathematical expectation

stand dev=sqrt(n .* p .* (1-p)) % standard deviation

m1=sum(x.*d) % mathematical expectation

s1=sqrt(sum((x-m1).ˆ 2.*d)) % standard deviation

pause

% comparison with normal distribution

hold on

dn=d gauss(x,math exp,stand dev);

plot(x,dn,’g’);

The function d=d bin(n,p) is represented below.

function d=d bin(n,p);

% binomial distribution d=d bin(n,p)

% n : n=1,2,...; - the maximum value of x

% x=0,1,...,n

% p : 0 ¡ p ¡ 1

n1=prod(1:n);

d=zeros(1,n+1); % determination of the row vector d

for m=0:n

if m==0

m1=1;

else

m1=prod(1:m);

end

nm=n-m;

if nm==0

Computer exercises 401

nm1=1;

else

nm1=prod(1:nm);

end

d(m+1)=(n1/(m1*nm1))*pˆ m * (1-p)ˆ nm;

end

A.1.3.3 Cauchy distribution

The density function of the Cauchy distribution is given by the formula (1.51).

Graphic presentation of the Cauchy density function is realized by the script

d0cauch.m. The function d=d cauch(x) return the row vector d, which contains

the values of the Cauchy density function.

A.1.3.4 χ

-distribution

Graphic presentation of the density function for χ

-distribution (1.45) is realized

in the script d0chisq.m. The function d=d chisq(x,n) returns the row vector d,

which contains the values of the density function.

The script d0chisq.m is represented below.

% Script d0chisq.m

% Calculation of the density function for χ

distribution

% call d=d chisq(x,n)

% n : a number of degree of freedom

% x : an argument of the density function, a ¡ x ¡ b

% d : returns the values of the density

% function

clear

% INPUT

n=50

% to determine an interval on x axis for computing the density

% function

a=0; % (a ¡ x ¡ b)

b=100.0;

nx=500; % sample size

% END INPUT

del=(b-a)/(nx-1); % to determine the values of argument

x=a:del:b;

d=d chisq(x,n); % row vector

plot(x,d); % graphic presentation

comp prob=trapz(x,d) % computing of the total probability

math exp=n % mathematical expectation

stand dev=sqrt(2 .* n) % standard deviation

402 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING

m1=trapz(x,x.*d) % mathematical expectation

s1=sqrt(trapz(x,(x-m1).ˆ 2.*d)) % standard deviation

pause

% comparison with the normal distribution

hold on

dn=d gauss(x,math exp,stand dev);

plot(x,dn,’g’); % graphic presentation

% of the normal distribution

The function d=d chisq(x,n) is represented below.

function d=d chisq(x,n)

% χ

-distribution d=d chisq(x,n)

% n : a number of the degree of freedom

% x : a row vector contains the values of argument (0 ¡ x ¡ inf)

% d : a row vector contains the values of the density function

d=((x./2).ˆ ((n./2)-1).*exp(-x ./ 2))/(2.*gamma(n./2));

A.1.3.5 Exponential distribution

Graphic presentation of the exponential density function (1.49) is realized by the

script d0expon.m. The function d=d expon(x,l) returns row vector d, Which con-

tains the values of the exponential density function.

A.1.3.6 Fisher distribution

Graphic presentation of the Fisher density function (1.48) is realized by the script

d0fish.m. The function d=d fish(x,n1,n2) returns row vector d, Which contains

the values of the Fisher density function.

A.1.3.7 Univariate normal (Gaussian) distribution

Graphic presentation of the normal density function (1.41) is realized by the script

d0gauss.m. The function d=d gauss(x,xm,s) returns row vector d, Which contains

the values of the normal density function.

A.1.3.8 Geometrical distribution

Graphic presentation of the geometrical density function (1.40) is realized by the

script d0geom.m. The function d=d geom(n,p) returns row vector d, Which contains

the values of the normal density function.