Troyan V., Kiselev Y. Statistical Methods of Geophysical Data Processing
Подождите немного. Документ загружается.
Computer exercises 403
A.1.3.9 Laplace distribution
Graphic presentation of the Laplace density function (1.50) is realized by the script
d0lapl.m. The function d=d lapl(x,l) returns row vector d, Which contains the
values of the Laplace density function.
A.1.3.10 Logarithmic normal distribution
Graphic presentation of the logarithmic normal density function (1.52) is realized
by the script d0lnorm.m. Function d=d lnorm(x,xm,s) returns row vector d, which
contains the values of the logarithmic normal density function.
A.1.3.11 Poisson distribution
Graphic presentation of the Poisson probability function (1.39) is realized by the
script d0poiss.m. Function d=d poiss(n,a) returns row vector d, which contains
the values of the Poisson probability function.
A.1.3.12 Student distribution
Graphic presentation of the Student density function (1.46) is realized by the script
d0stu.m. Function d=d stu(x,n) returns row vector d, which contains the values
of the Student density function.
A.1.3.13 Uniform distribution
Graphic presentation of the uniform density function (1.44) is realized by the script
d0unif.m. Function d=d unif(x,al,bt) returns row vector d, which contains the
values of the uniform density function.
A.1.3.14 Two-dimensional normal distribution
Graphic presentation of the two-dimensional normal density function (see Sec. 1.8.2)
is realized by the script d0gauss2.m. Function d=d gauss2.m returns matrix d,
which contains the values of the two-dimensional normal density function.
A.1.3.15 Exercises
(1) To compare the binomial distribution with the normal density function for the
following values of the parameters n (n = 1, 5, 20) and p (p = 0.1, 0.5, 0.9).
(2) To compare the Poisson distribution with the normal density function for the
following values of the parameter a (a = 0.1, 0.5, 1.0, 5.0).
(3) To compare a density function of the χ
2
-distribution with the normal density
function for the following values of the parameter n (n = 2, 5, 10, 20).
(4) To compare a density function of the Student distribution with the normal
density function for the following values of the parameter n (n = 2, 5, 10, 20).
404 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
(5) To implement a graphic presentation of the two-dimensional normal density
function for different mean values m
1
and m
2
(m
1
= 0, 1, 10; m
2
= 0, 1, 10),
standard deviations σ
1
and σ
2
(σ
1
= 0, 1, 10; σ
2
= 0, 1, 10), correlation coeffi-
cient ρ (ρ = −0, 95; 0; 0, 95).
A.1.4 Computer simulation of random values
There are the algorithms for simulation of the random values (see, for example,
(Ermakov and Mikhailov, 1982)). The scripts rnd0ald.m and rnd0alc.m return the
random values which belong respectively to discrete and continuous distributions
considered in Sec. A.1.3.
For these both cases a generator for the uniform random distribution on interval
[0; 1] with the following density function
f(x) =
1, 0 ≤ x ≤ 1,
0, x < 0 if x > 0.
(A.2)
is used.
For the discrete case the general method for the simulation of random variable
ξ with the probability function p
m
= P (ξ = x
m
) is based on the relation
P
m−1
X
n=0
p
n
≤ α <
m
X
n=0
p
n
!
= p
m
, (A.3)
where a value α belongs to the uniform distribution (A.2).
For the continuous case the formula for simulation of random variable using
cumulative distribution function F (x) can be written as
ξ = F
−1
(α), (A.4)
where F
−1
(α) is the inverse (cumulative) distribution function, α belongs to the
uniform distribution from (A.2) and a value ξ belongs to the distribution with
the cumulative distribution function F (x) (in (Ermakov and Mikhailov, 1982) the
various algorithms of the simulation of random values are considered). Because of
the strongly monotone continuous function F (x), we can write the next equalities
P (F
−1
(α) < x) = P (α < F (x)) = F (x). (A.5)
The scripts rnd0gau.m and rnd0uni.m return random values belonging to Gaussian
and uniform distributions correspondingly. Script rnd0n2.m returns a sampling
belonging to the 2-D Gaussian distribution.
A.1.4.1 Kolmogorov criterion
Script gc0kolm.m allows to realize the Kolmogorov criterion of belonging of the
sample cumulative distribution function to one or other type of the continuous
distribution. The Kolmogorov criterion gives probability of the random event, that
Computer exercises 405
a difference between cumulative distribution function F
ξ
(x) and its sample image
F
∗
N
(x) is smaller than given quantity γ/
√
N. In explicit form this criterion reads as
lim
N→∞
P
sup
−∞<x<∞
|F
ξ
(x) − F
∗
N
(x)| < γ/
√
N
=
P
∞
k=−∞
(−1)
k
exp(−2k
2
γ
2
), (A.6)
where N is sample size. Criterion is realized for continuous distributions from
Sec. A.1.3.
A.1.4.2 Exercises
(1) Using scripts rnd0alc.m and rnd0ald.m from Sec. A.1.4 to compare mean val-
ues and standard deviation with their population values for different sample
sizes (f.e. 10 – 500) for the cases of chi-square and Student distributions.
(2) Using script rnd0n2.m to investigate an influence of the sample size (f.e. 100
– 5000) on the sample mean and sample covariance matrix (correlation and
covariance coefficients).
(3) With the use of script gc0kolm.m to calculate the probability for the Kol-
mogorov criterion for the different types of the continuous distributions (Gaus-
sian, χ
2
, Student). To pay the main attention to a behaviour of the sample
(cumulative) distribution function in dependence of sample size (f.e. 10 – 1000).
A.1.5 Confidence intervals
Script g0conf.m returns the confidence interval for the sample mean of the random
value belonging to the Gaussian distribution (see Sec. 1.8.14).
A.1.5.1 Exercises
(1) Using script g0conf.m to find the confidence intervals for the sample mean of
the random value belonging to the Gaussian distribution for different values of
the confidence probability β (f.e. β=0.2, 0.5, 0.9) and sample size (10 – 100).
A.1.6 Time series
A.1.6.1 Correlated and uncorrelated time series
Script p0covrp.m returns uncorrelated time series and Markov chain (see Sec. 1.10),
creates sample covariance matrix.
A.1.6.2 Autocorrelation and crosscorrelation functions
Script p0corf.m returns autocorrelation function or crosscorrelation function com-
puted for correlated and uncorrelated processes.
406 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
A.1.6.3 Exercises
(1) Using script p0covrp.m to compare the covariance and correlation matrices for
correlated and uncorrelated time series. To use a different length of the time
series (f.e. 100, 1000, 10000).
(2) Using script p0corf.m to compare autocorrelation functions and crosscorrela-
tion functions for correlating and uncorrelated time series.
A.2 Transforms
A.2.1 Fourier transform
The MATLAB functions fft and ifft are realized the discrete Fourier transform
and inverse discrete Fourier transform respectively. These transforms are deter-
mined by the next relations (see also Sec. 12.1)
y
k
=
N−1
X
r=0
x
r
exp
−i
2πkr
N
, k = 0, . . . , N − 1 (A.7)
(Fourier transform),
x
r
=
1
N
N−1
X
k=0
y
k
exp
i
2πkr
N
, r = 0, . . . , N − 1 (A.8)
(inverse Fourier transform).
A.2.2 Signals and their spectral characteristics
Script p0wlet.m creates a simple seismic trace: wavelet and Gaussian noise. The
spectral characteristics, autocorrelation function, power spectrum of the trace are
calculated (see Sec. 1.10). The signals created by function p wlet.m, that allows to
use the simple signals in time domain.
(1) Discrete δ-function:
x
r
=
1, if r = j,
0, if r 6= j.
(A.9)
(2) Like δ-function signal:
x
r
= exp
−π(r − t
p
)/t
c
)
2
, r = 0, . . . , N − 1, (A.10)
where t
c
determines a width of the signal, t
p
determines the shift of the signal
relatively to the zero sample r = 0.
(3) Ricker wavelet:
x
r
= exp
−(π ∗ (r − (t
c
+ t
p
))/t
c
)
2
×
1 − 2 (π(r − (t
c
+ t
p
))/t
c
)
2
,
r = 0, . . . , N − 1, (A.11)
Computer exercises 407
where t
c
determines a width of the signal, t
p
determines the shift of the signal
relatively to the zero sample r = 0. An example of the graphic output of the
script p0wlet.m is represented at Fig. A.1.
Fig. A.1 An example of the graphic output ofp0wlet.m. Ricker wavelet (a), its spectrum module
(b), spectrum phase (c), real part and imaginary part of the spectrum (d). The abscissa axes
contain the samples (see formulas (A.7), (A.8)).
A.2.2.1 Exercises
(1) Using script p0wlet.m to calculate the spectra of signals (Ricker wavelet, δ-
function) for different values of t
p
and different signal-to noise ratio.
A.2.3 Multifactor analysis
Script p0cater1.m is realized the multifactor analysis (see Sec. 12.9).
Let us consider an example of the application of the script p0cater1.m to pro-
cessing of the simple trace “signal with noise” (see Fig.A.2). The input data are:
np=500; % length of the trace (sample)
jm=20; % maximum number of principal
% components
unzer=2; % number of principal components
% (nonzero eigenvalues)
408 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
% for recovery
op=3; % sinal shape
tp=100; % the time shift of the wavelet
xm=0; % expectation of the noise
s=0.05; % st deviation of the noise
The covariance matrix and its eigenvalues (and eigenvectors) are presented at
Fig. A.3 and Fig. A.4 correspondingly. The diagrams of the pairs of the principal
components and a complete set of the principal components are presented at Fig. A.5
and Fig. A.6 correspondingly (see Example from 12.9). The result of recovery of the
record using a different number (1 – 4) of the principal components are presented
at Fig. A.7.
Fig. A.2 Initial record without noise (a) and with additive uncorrelated noise (b).
Fig. A.3 Correlation matrix M = 20.
Computer exercises 409
Fig. A.4 Proper numbers (a) and eigenvectors (b) of the correlation matrix.
A.2.3.1 Exercises
(1) Using script p0cater1.m to implement processing of the trace (wavelet with
Gaussian noise) for the cases of different values of the maximum principal com-
ponents (5 – 20), different number of the principal components using for the
recovery (2 – 10), different values of signal-to-noise ratio.
410 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
Fig. A.5 The scatterplots of pairs of the principal components: M and M − 1 (a); M − 1 and
M − 2 (b); M − 2 and M − 3 (c); M − 3 and M − 4 (d). The maximum number of principal
components is 20.
Fig. A.6 The complete set of principal components (M = 20).
Computer exercises 411
Fig. A.7 Recovery of the record from Fig. A.2(b). Recovery using the first principal component
(a), two (b), three (c) and four (d) principal components (see Fig. A.6).
A.2.4 Cepstral transform
Script p0cep1.m is realized the cepstral transform (see Sec. 12.6).
The cepstral transform in the script p0cep1.m is realized as follows:
d=fft(trace); % Fourier transform
% of the vector “trace”
dabs=abs(d); % computing of the absolute value
% of the spectrum
dang=unwrap(angle(d)); % calculation of the phase
dcep=fft(log(dabs)+i*dang); % Fourier transform
% of the logarithm of spectrum (cepstrum)
The inverse cepstrum transform is realized as following:
dback=ifft(dcep); % inverse Fourier
dexp=exp(dback); % transform
dback0=ifft(dexp); % inverse Fourier
% transform
In Fig. A.8 and Fig. A.9 two signals and their cepstral characteristics are rep-
resented respectively. The result of the convolution of the two signals is compared
with the inverse cepstral transform of the sum of their cepstrums (see Fig. A.10).
412 STATISTICAL METHODS OF GEOPHYSICAL DATA PROCESSING
Fig. A.8 The first signal (a) and its cepstral characteristics (b) – (d). The Ricker wavelet (a),
real (b), imaginary (c) parts and the module of its cepstrum (d).
Fig. A.9 The second signal and its cepstrum characteristics (see Fig. A.8).
A.2.4.1 Exercises
Using script p0cep1.m to compare the convolution in time domain with summation
in cepstral domain.