Tip:
Highlight text to annotate it
X
In the previous lecture, we discussed methods for simulating samples of random variables
on a computer. So, we talked about pseudorandom number generators, that means, deterministic
algorithms which would help you to simulate random numbers, which are distributed uniformly
in 0 to 1. And subsequently, based on methods of transformation or accept reject methods,
we are able to simulate samples of scalar or vector random variables, Gaussian or non-Gaussian
random variables; and we also discussed what happens if random variables are completely
specified and what happens if they are partially specified.
So, this is what we discussed in the previous lecture. Now, in the present lecture, we will
now consider how to extend these capabilities to simulate samples of random processes.
Now, we will go back to the discussion that we had earlier on Fourier representation of
a Gaussian random process. So, we start with a 0 mean, stationary, Gaussian random process,
and we define in terms of random variables a n and b n as shown here. Suppose X of t
is n equal 1 to infinity a n cos omega n t plus b n sin omega n t omega n is n omega
naught; we assume that, this a n and b n are random variables, and a n is normal with mean
0 and standard deviation sigma n, b n is normal with mean 0 and standard deviation sigma n;
a n, a k are the independent for n naught equal to k; there uncorrelated, there Gaussian,
therefore independent. Similarly, b n, b k are independent for n naught equal to k, and
a n, b k are independent for every n and k. So, now, first let us deduce what are the
properties of this random process. What will be the mean of this process? Expected value
of X of t is, you have to take expectation inside, expected value of a n is 0, b n is
0, therefore, expected value of X of t is 0.
How about its covariance? Expected value of X of t and X of t plus tau. So, X of t is
the first summation here, n equal to 1 a n cos omega n t plus b n sin omega n t; the
second term n equal to 1 to infinity a n cos omega n t plus tau plus b n sin omega n t
plus tau. Now, this becomes a double summation, and if we expand and multiply, and use the
fact that, this assumptions on mutual dependents of a n and b n, a k and a n b k, etcetera,
we can show that, this auto covariance is indeed a function of only the time difference;
and indeed, we get R xx of tau has n equal to 1 to infinity sigma n square cos omega
n tau. So, therefore, it follows that x of t is has
0 mean and covariance function which is function of tau; therefore, it is a white sense stationary
process, but since X of t is Gaussian, because a n b n are all Gaussian and we are doing
a linear transformation on Gaussian random variables; X of t is also Gaussian. Therefore,
X of t is also strong sense stationary.
Now, let us consider we start the argument from a slightly different direction. Now,
let us consider a power spectral density function, which is given by n equal to 1 s omega n delta
omega n delta omega minus omega n. Corresponding to this, I get a covariance function which
let be call it as R tilted R xx tau, which is 1 by 2 pi; and the corresponding transformation
following the relationship between covariance and psd, and we can show that, R xx tilde
is 1 by 2 pi n equal to 1 S omega n delta omega n cos omega n tau.
Now, let us compare this with the auto covariance of X of t, that we just now obtain. We had
n equal to 1 sigma n square cos omega n tau; and in this case, if I put S omega n delta
omega n by 2 pi as sigma n square, you can see that the two representations are identical.
So, that would mean, if this is the target power spectral density function of X of t,
we can discretize the power spectral density into some n intervals; and the area under
this interval I call it as S omega n delta omega n. And if I consider now for each of
these intervals corresponding to the two random variables a n b n and use this representation,
we can show that samples of X of t, as n becomes large, we will be having power spectral density
which is the continuous form of the specified power spectral density function. So, this
gives as a approach to simulate samples of Gaussian random process is stationary 0 mean,
which specified power spectral density function. So, you start with the power spectral density
function, discretize, and for each segment, you define two random variables which are
independent and have the variance given by this area under the psd, and use that in this
representation, and you will able to produce sample samples of X of t, which has power
spectral density as given here.
So, let us see an example, simulate samples of a 0 mean stationary Gaussian random process
with properties that, S of omega is an as a as selected for illustration an exponential
type of a Gaussian type of power spectral density function; this is almost similar to
a Gaussian probability density function, except that, it is known a psd function; it is area
under the curve need not be 1, it will be equal to the specified variance, and some
numerical values I have provided for the parameters: I and alpha, which appear here. Now, X of
t I will write it as, i equal to 1 to n a n cos omega n t plus b n sin omega n t. So,
I am actually discretizing the power spectral density function into, each dot is a point
where I have discretized; and the parameters used are, we are simulating this for about
5 seconds, and I am retaining 120 terms and omega naught is taken us 0.2513 radian per
second; the omega max is 120 points is about 30 radian per second, and delta t which will
be function of 1 by omega max, is a 0.0419 seconds.
So, I simulate 240 Gaussian random variables with their mean being 0 and variance obtained
from the target power spectral density function. You have to find area under those small intervals
and you have to assign them as standard variances to the those one of, each one of those 240
random variables.
So, if I do that, one of the sample that I got is shown here - blue line; this is a sample
and 1000 such samples were obtained, and the ensemble mean and stand deviation were computed.
So, you see here, the target mean is 0 and the black line that you see houring around
0 is the simulated mean and this green line a pink line is a target standard deviation,
and the green line that is fluctuating about this is the simulated standard deviation;
so, things look quite.
And the probability distribution function with 1000 samples; blue is a simulation and
red is a target. There is reasonable agreement; again, we can do a hypothesis test to actually
objectively access, whether that agreement is acceptable or not.
And this is the ensemble of the time history; all of them I have shown on the same well.
This is again characterless, but gives a visual impression of what we are simulating.
Another example, we discussed this Kanai-Tajimi power spectral density models, when we discussed
a models for earthquake ground accelerations. Suppose I want to simulate now samples of
earthquake ground accelerations, whose power spectral density function is given by the
Kanai-Tajimi power spectral density function, we some of these parameters, I equal to 1,
omega g is the ground natural frequency, eta g is the ground damping. So, the idea here
is that, the soil layer is underlying a overlying a bedrock is modeled as single degree freedom
system, with damping eta g and natural frequency omega g, and it is subjected to a white noise
acceleration at the bedrock level. And the ground surface acceleration is modeled as
output of this single degree freedom system in steady state to the white noise excitation
applied at the bedrock level.
So, the target power spectral density function is the blue line axis, y axis is shown on
large scale and x axis is on linear scale here. So, after doing I think 5000 simulation,
we are able to estimate the power spectral density from the data and it is compared with
the target, and again, we see good agreement between the 2. So, in a future lecture, I
will describe how to actually estimate power spectral density from data, and what are the
sampling distribution for power spectral density functions and how to test hypothesis, etcetera
will do that later, but right now will be a happy with a kind of visual comparisons,
and this that level, the results are satisfactory.
This is a display of the a few samples of the time histories. Again, to give a visual
idea of the how the samples look, we are simulating samples of stationary random processes for
a chosen time length. This is the an example where... this another example where the covariance
function is a harmonically decaying - exponentially decaying harmonic - with certain parameter
displayed here and target density function is Gaussian. And what is shown here is the
results on probability distribution function; blue is the simulation, red is the normal,
which is the target just has on a illustration.
This is for the same example; this is auto covariance function. Its Fourier transform
will give us a power spectral density function, and one sided power spectral density function
is shown here; blue is the target and red is the simulations obtained using 500 samples.
So, there will be fluctuations - sampling fluctuations - that arise, we you need to
understand that.
This is a result on probability distribution function; again, the comparisons are satisfied;
some samples of X of t. Now, so simulation of a scalar Gaussian random process with,
which is stationary with 0 mean and given auto covariance are equivalently the power
spectral density function, seems reasonably straight forward. Non-stationarity can be
introduce by multiplying stationary processes by suitable envelops, that is one of the common
strategy is used; we will discuss that later. But right now, we will move on to the problem
of simulating samples of scalar non-Gaussian random processes. The description of the random
process here is a limited to the power spectral density function or the covariance function
- auto covariance function - and first order probability density function. So, here, we
are considering, let X of t be a random process whose first order pdf and auto covariance
functions are available; no further information about the process is available is taken to
be available. X of t need not be stationary. So, the problem is how to simulate samples
of X of t. So, what we define? We follow the Nataf's
of transformation method that I discussed in the previous lecture. So, first step is,
we remove the mean and divide by standard deviation, so that Y of t has 0 mean and unit
standard deviation. Now, you introduce a new random process Z of t through the transformation,
phi of Z of t is equal to P y of Y of t, where phi of this argument is probability distribution
function of a normal random variable with 0 mean and unit standard deviation. Z of t
is a 0 mean Gaussian random process with an unknown covariance matrix or the covariance
function here.
So, we follow the argument, we have shown in the previous lecture that, as for the first
order probability distribution function is concerned, the properties are preserved; so,
we not have to go through the step again. So, samples of Y of t can be simulated following
this rule, one samples of Z of t are obtained. Now, I need the covariance function of Z of
t; so, what I know is a covariance of Y of t. So, I will express this in terms of the
covariance of Z of t, this rho star t 1 comma t 2 is not known; what is known is here, this
is known, this is unknown. So, as the prelude to the implementation of
this method, we have to follow that approach, that I described in the previous lecture.
We know, although with this is not known, we know that this lies between minus 1 and
plus 1; the left hand side also this rho xx is also takes values in minus 1 to plus 1,
this is known, this is unknown. So, what we do is, we find out the left hand side for
various specified values of rho star from minus 1 to plus 1 and compute the corresponding
rho; and among these computed values of rho, there will be one value which is the target
value. So, we will be able to find out the corresponding rho.
Now, this step can be implemented; so, this step, the steps involve in simulation would
be solved for rho star t 1 comma t 2 and simulate z of t, and then simulate Y of t using this
rule, and then go back and simulate X of t by using this rule. So, the computational
steps involved are quite tedious, but in principle, it is dual.
So, things work out become bit simple, if a process is stationary. So, indeed I start
with the case of a stationary random process, whose covariance is sigma square exponential
minus alpha tau and parameters of that is given; and I want the density function to
be uniformly distributed between minus 0.5 plus 0.5. So, first we have to find out rho
star, and then do this transformation - Nataf's transformation - and this some results are
displayed here; red is a target and blue is a simulation, and the agreement is quite satisfaction.
Now, this is a target power spectral density function, which is actually the Fourier transform
of this, which is something like sigma square alpha by alpha square omega square, something
like that. So, that is displayed here; the red and the blue, blue is the target and red
is the simulated estimated one; they show quite a good agreement and we are unable to
distinguish the two curves on these plot, and this is obtained with 500 samples.
These are few samples of X of t, and notice that, that is uniformly distributed between
minus 0.5 to plus 0.5 and it has 0 mean. And if you compute the power spectral density
function or auto covariance function, it will match with the target values - target for
values - that is specified in the statement of the problem.
Same auto covariance function but will slightly different numerical values, but now the density
function is not uniform, it is Rayleigh. So, mean would not be 0, the samples will have
non-zero mean; and the parameters are specified here, sigma that appears in the definition
P x of x is 2. How do you a simulate samples of this process? So, again, we follow this
Nataf's transformation method and go through that. Again, let me emphasize this description
of the random process is not a complete specification of the random process; we are only telling
the marginal density function and a covariance function. If the process is Gaussian, this
is adequate; the first order covariance will be sufficient; for characterize second order
properties and the first order property would give you the mean, second order is covariance;
that is enough to completely specified Gaussian random process.
But we are talking about non-Gaussian random processes. So, this is not a complete specification.
So, we are able to simulate samples, whose properties match with what has been specified;
and other properties of this process, we really do not know what it correspond to it; it is
left to the implicit modal that is contained in the Nataf's transformation.
Now, this is the power spectral density function further Rayleigh random process. And the power
spectral density function is shown on a log-log scale to make the comparisons clearer, and
again, the agreement between target which is blue and simulation which is red is quite
ok.
This is samples of a Rayleigh random variable and the psd that we mentioned that random
process. And again, we can see that mean of these samples are non-zero, it is oscillating
apart some non-zero value, which is what we should be expecting for this model.
Now, scalar random processes we are able to simulate Gaussian, we have to discuss target
power spectral density or target covariance, we are able to match. And for scalar non Gaussian,
the target first order PDF and the specified PSD are auto covariance we are able to match.
Now, how about vector Gaussian random processes? So, to illustrate that, let us consider two
random processes X of t and Y of t; this is specified to be jointly stationary and they
have 0 mean. So, the covariance - auto covariance and cross covariance functions - completely
specify the random process. So, the auto covariance of X of t is expected
value of X of t into X of t plus 2, which is R xx tau. Similarly, R yy of tau is defined
like this; R xy of tau is the cross covariance function between X of t and Y of t which is
expected value of X of t and Y of t plus tau. Now, let us look at R xy of tau, this is X
of t in to Y of t plus tau. So, this is actually we can this is equal to Y of t plus tau in
to X of t; therefore, this should be equal to Y of t plus tau in to X of t and this is
nothing but R yx of minus tau. So, R xy of tau is R yx of minus tau; so,
it is not a symmetric function, but there exist certain skew symmetric of this kind.
Now, R xy of tau is X of t into Y of t plus tau, but this is not equal to Y of t in to
X of t plus tau; therefore, R xy of tau is not same as R y of x of tau. So, the R of
tau is specified in terms of these functions; it is not symmetric, but there exist certain
relationship between the two.
Now, how about the corresponding descriptions in frequency domain? So, the Fourier transform
pair of R xx and S xx are displayed here; the R xx expressed in terms of S xx and S
xx this expressed in terms of R xx. So, similarly R xy is shown here. Now, R yx of tau, we can
write in this form and defined S yx of omega; and since R xy of tau is not symmetric, they
corresponding power spectral density function will be complex valued, and it has real part
and an imaginary part. So, we can show that, since R xy of tau is R yx of minus tau, it
follows that S xy of omega is same as conjugate of S yx of omega. So, these relations exist;
there not four independent functions, but they are not symmetric either. So, there is
a interdependence between the four quantities in frequency and time domains.
Let us examine this cross covariance and cross power spectral density function slighter slightly
greater detail. So, S xy of omega is gamma xy of omega plus i delta xy of omega. Now,
substitute this into the definition R xy of tau and I get this gamma plus i delta, and
for e raise to minus i omega t, I will write it is as cos omega t minus i sin omega t;
I can separate real and imaginary parts. Now, although S xy of omega is complex valued,
R xy of tau is real valued. So, that would mean the imaginary part of this must be equal
to 0; therefore, this second integral must be 0. From this, it follows that gamma xy
of omega is an even function and delta xy of omega is an odd function, that ensures
that R xy of tau is real valued.
So, this is the summary of the relationship between R xy and gamma and delta. And we can
use that property and write R xy of tau as 0 to infinity 1 by pi now, not minus infinity
to plus infinity. See, gamma xy is even, cos is even, the so the product is even; this
is odd, sin is odd, so odd into odd is even. So, both this integrant here is even; therefore,
I can write it as 2 into 1 by pi 0 to infinity.
Now, the various definitions of power spectral density function in terms of an expectation
of a truncated Fourier transform; this is well known; we have already studied this.
And based on this, we can write the matrix of power spectral density functions, and S
xy of omega and S yx of omega are not the same, but nevertheless, they are related through
this relation, that S xy of omega is S star yx of omega.
Now, let us do start with Fourier representation for X of t and Y of t. I will introduce for
X of t, two families of random variables a n and b n; for Y of t, two families c k and
d k. Now, a n and b n, see, the X of t is stationary, Y of t is stationary, we would
as a scalar random process. So, a n whatever conditions we impose on a n and b n, whenever
representing X of t will now must continue to apply; that would mean a n is normal 0
mean and sigma xn, b n is similarly 0 mean sigma xn standard deviation, a n a k are independent
for n not equal to k, b n and b k are independent for n not equal to k, and a n b k are independent
for all n and k.
So, X of t as for as a scalar random process is concerned; it has 0 mean and its auto covariance
is given by this. So, this is stationary random process, this we have seen. The same logic
can be used on Y of t, and we again we can show that, you know, c n is normal 0 mean,
sigma y n is standard deviation, and so on and so forth. And expected value of Y of t
is 0, covariance of Y of t is n equal to 1 sigma y n square cos omega n tau. So, this
sigma xn square and sigma y n square are to be obtained from discretize version of auto
power spectral density function of X and Y respectively, so that we have already seen.
Now, what is more what is now remains to be addressed is the properties cross covariance
properties and cross power spectral density function properties and capsulated in those
functions. So, how do we deal with that?
So, let us consider the expectation of X of t into Y of t plus tau. So, for X of t, I
have this representation; for Y of t, I have this representation. I will use this and write
it as the first term is this first summation, and the second term is the second summation,
where t is replaced by t plus tau. So, I run through this calculations, the product of
summations become double summation and I will now use the properties that we have encapsulated
for a n and b n.
What we will do now is, I have not defined the properties between a n and c k how they
are related, a n and d k how they are related, b n and c k how they are related, b n and
d k how they are related, that we have not specified. So, what do will do is, we will
assume that, a n is independent of c k and d k, if n is not equal to k; similarly, b
n is independent of c k and d k, if n not equal to k; that we will be it, but when n
equal to k, I will impose certain restrictions on expectations of a n, c n, a n d n, b n
c n and b n d n, so that I will be able to simulate the properties of contain properties
contained in the cross covariance functions. So, to do that, what I do here? First of all,
I want that, the properties that X and Y are jointly stationary must be honored. So, if
that has to be honored, the expected value of X of t in to Y of t plus tau must be function
of tau alone. So, for that have to happen, if I take now expected value of a n c k is
sigma a c n, delta nk, this is a chronicle delta, which is equal to 1, if n equal to
k; otherwise, it is 0. And similarly, if I place restrictions on a n d k, expected value
of a n d k, expected value of b n c k, and so on and so forth, I can show that the expression
for R xy can be reduced; the double summation collapses into a single summation.
Now, I will make further assumption, that sigma acn is minus sigma bdn and sigma adn
is minus sigma bcn. If I do that, I can show that R xy of tau now becomes a function of
tau alone. Now, I can select sigma acn and sigma adn to make sure that, this representation
actually corresponds to the target value of cross covariance functions.
How do we do that? Now, we will consider now R xx, R yy, R xy; this is the 1, 2, 3 are
the consequence of are the assume properties of a n, b n, c n, d n, for n from 1 to capital
N. S xx omega if I consider S xx omega given in this form, we already seen that the corresponding
auto covariance function matches with, can be made to match with function of this form
by selecting sigma x n square to be equal to S xx omega n delta omega n divided by 2
pi.
So, that means, this is the target and this is what happens if I use the power spectral
density that I mention just now. And these two can be made to agree among themselves,
if sigma x n square is taken to be this quantity. So, that R xx tau becomes R tilde of xx of
tau. Now, similarly, S y of omega I define in a series like this, the power spectral
density functions; and if I select sigma y n square to be given by this, I can show that
the Fourier transform of this, which is R tilde yy of tau will match with the target
auto covariance of Y; this is similar exactly similar to what we did for X of t.
How about R xy? R xy of tau is given by this, and by separating real and imaginary parts,
and honoring the fact that R xy of tau is indeed real valued, I get this functions.
I have already shown that the imaginary part is 0, and what is the implication of that
on properties of gamma and delta. Now, if I now consider gamma to be a sequence like
this and delta to be a sequence like this, and select, of course, you corresponding to
this sequence, you can take a Fourier transform; and that is quite easy, because we have direct
delta functions the integration is very straight forward. If I do that, the expression for
the cross covariance function can be shown to be given by this form; it is again of the
form some constant of a n multiplied by cos omega n tau constant function of n multiplied
by sin omega n tau. Compare this with R xy of tau, and that is suggest that, if I select
sigma acn to be this and sigma adn to be this, the cross covariance are tilde xy and R xy
of tau match.
So, I have now the recipe to simulate the vector Gaussian random process. So, summary,
I use this Fourier representation for X of t, Fourier representation for Y of t; this
sigma x n that appears as standard deviation of a n is obtain from psd of x, sigma y n
and that appear for c n and d n is obtained from power spectral density of Y of t, and
sigma acn and sigma adn that appear in the description of properties between a n c n,
a n d n, b n c n and b n d n is obtained from gamma and delta of your cross power spectral
density function. So, the out of power spectral density function of x and y, the gamma and
delta functions for correspond to the cross power spectral density between x and y are
given. So, from that, I can always by discretizing
them, I can always get sigma x n and sigma y n, sigma can, sigma adn and we are ready
to use this series. Once we are able to simulate a n, b n, c n, d n according to this prescription,
we will be able to simulate x n y n.
So, summary is a n, b n, c n, d n are four-dimensional normal random variables with 0 mean and covariance
given by this. So, to implement the method, we should be able to simulate samples of a
n, b n, c n, d n, according to this prescription. So, here, again you need to do the calculation
that I mentioned, namely, you have to find the eigenvalues and eigenvectors of this,
and do the transformation, and first in the standard normal space and then work backwards
and get a n, b n, c n, d n, where all these parameters appearing here are expressed in
terms of known properties of X of t and Y of t.
So, after that, we get expression for R xx, R yy, R xy, which is quite satisfactory. Some
more notes, expected value of X of t into Y of t plus tau is R xy of tau, and expected
value of Y of t plus tau in to X of t is R yx of minus tau, which is same as R xy of
tau. Now, if you consider expected value of Y of t minus tau in to X of t, which is R
yx of tau, this is same as R xy of minus tau. So, all these you can verify. Now, in our
calculation, we got sigma x y of tau a sigma acn cos omega n tau sigma adn sin omega n
tau sum from 1 to n, which is not equal to R yx of tau which is R xy of minus tau, because
sin omega n of minus tau is minus of sigma adn sin omega n tau. So, things are ok, these
are some simple checks on basic properties, that we should expect from such representations.
Now, I leave it as an exercise; this is the fairly... this can be a fairly long exercise
which can be a accomplished, if you write a computer program. So, as an illustration,
we will consider as simulation of spatially varying earthquake ground acceleration. Consider
X of t and Y of t to be two random processes, representing ground accelerations in the horizontal
direction at two stations separated by a distance d x y; that means, the same event of an earthquake,
I take two points on the earth surface and consider the ground acceleration in a given
direction; and at one point, I call it as X of t, and another point, it is Y of t.
We can take that X of t and Y of t to be jointly stationary, Gaussian and 0 mean random processes.
The auto psd functions of X of t and Y of t may be taken to be of the form, S of omega
is I in to H 1 of omega whole square and H 2 of omega whole square; this is typical Kanai-Tajimi
power spectral density function model, where H 1 is given by the transfer function of the
soil layer and H 2 of omega is an artificial filter function, which eliminates certain
anomalies in the low frequency behavior of samples of displacements, that is needed to
make the model work well at low frequencies.
So, the auto psd's are specified; this is the well-known Kanai-Tajimi power spectral
density function model. And we introduce a coherency function, that is, S xy by S x square
root S xx S yy to be given by this; this is a complex valued coherency function, which
depends on the distance between the two stations, and the wave velocity and frequency, and it
is a complex valued function. So, from this, you can actually get multiply this by the
2 power spectral density functions, you get the cross psd function; and you can separate
them into real and imaginary parts, you will get your gamma and delta; so, you are ready
to launch your simulation.
So, problem on hand is to develop a computer code to simulate samples of X of t and Y of
t. And to facilitate that exercise, some numerical values I have suggested here; ground natural
frequencies 15.6 radian per second, and so on and so forth. So, the exercise consist
of simulating, may be say 5000 samples of X of t and Y of t, and then follow it by...
From that ensemble of time histories, estimate the power spectral density function matrix
and compare it with the target psd matrix. So, this is a quiet an involved exercise.
This second part of this exercise, namely, estimating psd matrix from observe data of
vector realizations, is something that we need to consider later. At this stage, we
do not have I am not discuss this yet, but it is a part of this exercise.
Now, I have talked about pairs of random processes, we can generalize it to say three random processes
or n random process. So, this method seems to work out quite well for this. So, we use
three Fourier representations, and I have now auto psd for x, auto psd for y, auto psd
for z, cross psd between x and y, x and z, and y and z. So, I have to now each one is
a stationary random process; so, represent finding properties of a n and b n is purely
based on property of psd's of x, y and z, for say, a n b n, c k d k, e m f m, but to
find the cross properties between a n c n, a n e n, etcetera. You have to bank on the
cross power spectral density functions; so, you can specify all this properties, as for
the formulation that have just now outline; and in principle, it should be possible to
simulate samples of n dimensional Gaussian vector random processes.
I have detail the various properties that we need to use and I leave that exercise for
you to verify. How about multi parameter random processes? For example, if you are modeling
wind velocity on a say a chimney, so the wind velocity varies in time, but also varies the
space. So, I have a random process which is a function
of x n t. so, in a three-dimensional situation, this x can be x 1, x 2, x y z and t, like
ocean waves, and so on and so forth. But here, I am considering a single random processes
function of two parameters. This is a scalar situation; so, we can a vector situations
also here; instead of just velocity, you can model some other parameter like pressure or
something like that. You have two quantities, but we will consider scalar case first. We
will assume that, this random field has 0 mean and it is homogenous, homogenous of homogeneity
of random field is synonyms are analogous to stationarity of a random process, that
would mean expectation of f of x comma t into f of x plus i plus t plus tau is a function
of x i and tau, that is the stationarity property. And corresponding to the two-dimensional auto
covariance function, I can define a two- dimensional power spectral density function. One of the
parameter will correspond to time, other one to space. This is the wave length, this is
the frequency in radian per second, this will be wavelength per unit length. So, I can define
the power spectral density function. Now, I can use a Fourier representation involving,
A nk, B nk, C nk, D nk and products of cos and sin functions, evolving in sin and evolving
in time and space. Now, these are 0 mean Gaussian random variables; I will adjust properties
of these random variables to suit the properties of this specified power spectral density function.
So, the problem on hand is, to simulate samples of f of x comma t, whose power spectral density
function in lambda and omega space is given; it is and the process is given to be Gaussian.
So, you can set up the same logic; we can make all these things for n naught equal to
k etcetera to be independent. And when for certain cases, when n equal to R n, k equal
to S, etcetera, we will adjust those parameters and take them to correspond to the value of
the discretize version of the power spectral density function matrix, and we can show that
they auto covariance from such a model. And from a model which is purely based on discretize
version of power spectral density function would match, if this A nk B nk are selected
in a certain particular manner and you can simulate samples.
Now, then alternate approach for simulation of non-Gaussian random processes with prescribe
power spectral density function. This is based on Markov process theory. So, let us consider
this problem; let X of t be a stationary Gaussian random process defined on interval x l to
x r, let P x of x be the first order probability density function, and let the psd function
for purpose of illustration, I am taking it to be of this form. Now, I want to simulate
the samples of X of t; now, what I do is, I consider a stochastic differential equation,
dX of t is minus alpha Xdt plus D of x dB t.
Now, this equation is such that, the two time response of X of t can be evaluated and we
can in fact compute the auto covariance of the response process. The parameter alpha
and this parameter D of X are unknowns; although I have selected alpha to be same as this,
it indeed transverse, those two will be the same, but the problem statement is as follows.
Suppose this alpha is say alpha star, so the drift and diffusion coefficients are unknown.
For this problem, the steady state solution to the governing of Fokker Planck equation
is obtainable; first order equation I showed already, that you can solve the Fokker Planck
equation in the steady state. Now, what I do is, I demand that if I where
to obtain the steady state solution of this equation on first order probability density
function, you should match with this P x of x. And similarly if I compute the power spectral
density function in steady state, it should match with this. So, what is known here, are
the solutions of this problem. What is unknown are this alpha star and D of X; so, these
are inverse problem. Now, what we do is, we demand that these solution match with the
target psd and pdf; that means, determine the drift and diffusion coefficients, so that
this becomes possible, that means, we complying with the target psd and pdf, and it becomes
possible.
So, how do we do that? So, we can start by multiplying the first of the question by X
of t minus tau and take the ensemble average. And from this, we get the equation for covariance
function in steady state, and we can show that, the covariance function is of the form
a exponential minus alpha mod tau. And corresponding to that, I have this power spectral density
function; that means, this alpha I should I have written as alpha star, but so imagine
this alpha star, this is alpha star; as for as psd function is concerned, select alpha
equal to alpha star. Then, our target is on psd in met; the target on first order probability
density function is yet to be satisfied. Now, what you should notice is, that the diffusion
coefficient D of X has no influence on psd of X of t in this equation.
Now, let us consider now the governing Fokker Planck equation. So, dX of t is minus alpha
Xdt plus D of X dB t, and associated with that, this is the governing Fokker Planck
equation. In the steady state, dou p by dou t would be 0 and I get this reduced equations.
So, one of the solutions for this would be of this form, and indeed, we can solve this;
this is the first order ordinary linear differential equation and we can solve this. Here, what
is not known? Alpha is known, because I already I solved for that; P of x is given, right;
D of x is unknown. So, this has to be viewed as an equation for D of x; that means, for
what value of D of x, would P of x would be the solution to this problem? Before so, it
is an inverse argument; so, if I do that and solve for D square of x, I get D square of
x to be this. That means, if I consider now a stochastic differential equation, whose
drift is alpha and diffusion coefficient is D square of x, square root of this is given
by this, and simulates samples of X of t according to this rule, the samples in steady state
would have this power spectral density and the target P of x as a probability distribution
function. So, that means, where Taylor the SDE is Taylor
made, to produce the requisite power spectral density function and first order probability
density function; it is an inverse approach. So, he has an exercise, you can consider producing
a samples of Rayleigh random process with this power spectral density function or uniform
distribution, I leave it as an exercise; I will not be able to, I will not discuss the
details.
Now, I would like to briefly introduce a new topic, what is known as variance reduction.
So, before that, we can summarize what all we have done. We have develop methods for
simulation of random variables, scalar, vector, Gaussian, non-Gaussian, completely specified,
partially specified. Similarly, methods for simulation of samples of random processes,
scalar Gaussian, scalar vector, scalar non-Gaussian, vector non-Gaussian with partial specifications,
and also discuss briefly a method based on Markov processes characteristics of the solutions.
So, as for as simulation tools are concerned, we same to be having the requisite a good
deal of tools already available with us.
Now, the problem of dynamic response characterization has to be addressed now. So, what we have
done is, we have simulated samples of actions. Now, we have to run them through the system
dynamics and get the samples of the response, and then be able to solve them and get properties
of response characteristics. So, one of the problem that arises there is a is known as
variance reduction, and I will briefly explain this in the broader context of Monte Carlo
simulations. And I will consider the problem of evaluation of a multi-dimensional integral
given by, integral P x of x dx over a region g of X less than or equal to 0.
So, this is the problem that arises in structural reliability analysis. So, the details of this
statement of the problem and its relation to problems structural reliability is not
of primary concern here. What is of concern here is the computational details of evaluation
of this integral using Monte Carlo simulations. So, I can, now, the region over reach integration
is being done appears as a limits of this integral, I can first make the limits and
minus limit to plus infinity by introducing an indicator function on g of x; this I of
g of X is 1, whenever g of X is less than 0; otherwise, it is 0. So, I can write it
as I of g of X P x of x dx. Now, therefore, P f can be written as expected
value of indicator of g of x; therefore, this integral formally I am writing it as an expectation.
Now, my problem is to estimate this through samples of using samples of x. So, let theta
be the estimator, i running from 1 to n a i I g of X i. Now, the mean value of theta
is given by, i running to 1 a i expected value of this, and this is known to be P f and this
is i equal to 1 to n a i.
Now, if this summation is equal to 1, we get theta to be unbiased; variance of theta I
can do a simple calculation, and I can show that, variance of theta is actually given
by this equation, I equal to 1 to n a i square P F into 1 minus P F.
Now, a i is there still unselected. So, what we do? Select a I, i equal to 1 to n, such
that, the variance of theta is minimize subject to this constraint; that means, moment this
equation this constraint is satisfied, theta become unbiased estimator; and moment this
optimization criteria is met, we get an unbiased estimator with minimum variance.
So, we can solve for a i and we have gone through this step in the discussion on estimator
for the mean. And we can show that, a k is 1 by k, these are the optimal values of the
parameters a 1, a 2, a 3, a n and associated variance is given by square root P F 1 minus
P F by n.
So, this result is quite similar to what we did earlier, so in characterizing the estimator
for the mean. Now, let us look at this slightly more carefully. Now, look at the standard
deviation which is square root of the variance of the sampling distribution, we get it in
this form; and the coefficient of variation can be defined as sigma by m, that is, 1 by
P F into this and a slight rearrangement will show as, that this coefficient of variation
can be given as, 1 divided by square root P F n, when P F n is small, which is most
of the cases the kind of problem that we are studying.
Suppose the answer is in the range of 10 to the power of minus 5 and you are looking for
coefficient of variation of 10 percent; that is, accuracy with which you want compute the
value of the integral, then the number of samples needed becomes 10 to the power of
7; that means, you need 10 to the power of 7 samples to evaluate P F, where coefficient
of variation is about 10 percent. This coefficient of variation can be viewed as some kind of
an error. On the other hand, if coefficient of variation becomes 0.01and probability of
this P F, which is probability of failure in structural reliability language is of the
order of 10 to the power of minus 5; the number of samples needed becomes 10 to the power
of 9.
Now, we can make few observations. Now, these variance of estimator, which is P F in to
one minus P F by n is independent of size of the basic random variable x; x can be one
dimensional, ten-dimensional, hundred-dimensional, this expression is independent of that; this
n is samples not a dimensions of x. If this variance is large, the utility of estimator
becomes questionable, because our population parameter is deterministic; and we are using
a random variable to approximate an deterministic quantity, its variance should be as slow as
low as possible; so, larger variance means the answer is less acceptable.
So, the question of reducing the variance is very important. So, how do you reduce the
variance? You look at the expression for the variance, the n is in the denominator here;
so, it appears that, in order to reduce the variance of the estimator, we need to increase
the sample size n. But, is it the only approach? We can raise this question. Can you reduce
the variance of the estimator without increasing n? So, this problem is known as problem of
variance reduction. Variance can be increased by increasing sample size, but are there any
other ways of reducing the variance without increasing sample size?
So, for a given value of n, instead of using this as my estimator, can I modify this estimator
and reduce the variants? The sample size is fixed, ok. So, this question is of a considerable
significants, when we analyze the response, because if you are applying simulation methods
for practical structures, if one run of computer simulation of a dynamical behavior of the
system takes about, say 5 minutes of C P U, and you are evaluating probability of failure
of the order of 10 to the power of minus 5 with coefficient of variation of, say 0.01
as I was mentioning, you need to run the computer code 10 to the power of 9 times. And if 5
minutes in to 10 to the power of 9, it is an hopelessly large amount of computational
time and that approach is unlikely to work. So, we need to develop some intelligent methods
to reduce the variance and that issue will consider subsequently. So, at this point,
I just want to state, what is the problem of variance reduction. So, at this stage,
we will close this lecture. In the next lecture, we will consider questions on simulation of
dynamical behavior of systems; we have now completed the problem of simulating samples
of random variables and random processes, and this description will be useful in characterizing
the inputs and initial conditions to the say the dynamical system. And now, the question
of, how does these ensemble of random variables and random processes transmit through the
dynamics of the system and produce the response quantities of interest? So, we will address
this question in the forthcoming lectures.