Reading for Today's Lecture: ?
Goals of Today's Lecture:
Central Limit Theorems
Slutsky's Theorem: If Xn converges in distribution to Xand Yn converges in distribution (or in probability) to c, a
constant, then Xn+Yn converges in distribution to X+c. More
generally, if f(x,y) is continuous then
.
Warning: the hypothesis that the limit of Yn be constant is essential.
Definition: We say Yn converges to Y in probability if
The fact is that for Y constant convergence in distribution and in probability are the same. In general convergence in probability implies convergence in distribution. Both of these are weaker than almost sure convergence:
Definition: We say Yn converges to Y almost surely if
The delta method:
Suppose a sequence Yn of random variables converges to some y a
constant and that if we define
Xn = an(Yn-y) then Xn converges in
distribution to some random variable X. Suppose that f is a
differentiable function on the range of Yn.
Then
an(f(Yn)-f(y)) converges in distribution to
.
If
Xn is in Rp and f maps Rp to Rq then
is the
matrix of first derivatives of components of f.
Example: Suppose
are a sample from a population with
mean
,
variance
,
and third and fourth central moments
and
.
Then
We take Yn to be the vector with components
.
Then Yn converges to
.
Take
an = n1/2. Then
Remark: In this sort of problem it is best to learn to recognize that the
sample variance is unaffected by subtracting
from each X. Thus
there is no loss in assuming
which simplifies
and a.
Special case: if the observations are
then
and
.
Our calculation has
Monte Carlo
The last method of distribution theory that I will review is that of Monte Carlo
simulation. Suppose you have some random variables
whose joint distribution is specified and a statistic
whose distribution you want to know. To compute
something like P(T > t) for some specific value of t we appeal to the
limiting relative frequency interpretation of probability: P(T>t) is the
limit of the proportion of trials in a long sequence of trials in which Toccurs. We use a (pseudo) random number generator to generate a sample
and then calculate the statistic getting T1. Then we
generate a new sample (independently of our first, say) and calculate
T2. We repeat this a large number of times say N and just count up
how many of the Tk are larger than t. If there are M such Tkwe estimate that
P(T>t) =M/N.
The quantity M has a Binomial( N,p=P(T>t)) distribution. The standard error of M/N is then p(1-p)/N which is estimated by M(N-M)/N3. This permits us to guess the accuracy of our study.
Notice that the standard deviation of M/N is
so
that to improve the accuracy by a factor of 2 requires 4 times as many
samples. This makes Monte Carlo a relatively time consuming method of
calculation. There are a number of tricks to make the method more accurate
(though they only change the constant of proportionality - the SE is still
inversely proportional to the square root of the sample size).
Most computer languages have a facility for generating pseudo uniform random numbers, that is, variables U which have (approximately of course) a Uniform[0,1] distribution. Other distributions are generated by transformation:
Exponential:
has an exponential distribution:
Random uniforms generated on the computer sometimes have only 6 or 7 digits
or so of detail. This can make the tail of your distribution grainy.
If U were actually a multiple of 10-6 for instance then the largest
possible value of X is .
This problem can be ameliorated by
the following algorithm:
Normal: In general if F is a continuous cdf and U is
Uniform[0,1] then
Y=F-1(U) has cdf F because
This is almost the technique we used above for the exponential distribution.
For the normal distribution
(
is a common notation
for the standard normal cdf) there is no closed form for F-1.
You could use a numerical algorithm to compute F-1 or you
could use the following Box Müller trick. Generate U1,U2
two independent Uniform[0,1] variables. Define
and
.
Then you can check using the
change of variables formula that Y1 and Y2 are independent
N(0,1) variables.
Acceptance Rejection
If you can't easily calculate F-1 but you know f you can
try the acceptance rejection method. Find a density g and a constant
c such that
for each x and G-1 is computable
or you otherwise know how to generate observations
independently from g. Generate W1. Compute
.
Generate a uniform[0,1] random variable U1 independent of all
the Ws and let Y=W1 if
.
Otherwise get a new W and a
new U and repeat until you find a
.
You
make Y be the last W you generated. This Y has density f.
Markov Monte Carlo
In the last 10 years the following tactic has become popular, particularly
for generating multivariate observations. If
is an
(ergodic) Markov chain with stationary transitions and the stationary
initial distribution of W has density f then you can get random
variables which have the marginal density f by starting off the Markov
chain and letting it run for a long time. The marginal distributions
of the Wi converge to f. So you can estimate things like
by computing the fraction of the Wi which land in
A.
There are now many versions of this technique including Gibbs Sampling and the Metropolis-Hastings algorithm. (The technique was invented in the 1950s by physicists: Metropolis et al. One of the authors of the paper was Edward Teller ``father of the hydrogen bomb''.)
Importance Sampling
If you want to compute
Variance reduction
Consider the problem of estimating the distribution of the sample mean
for a Cauchy random variable.
The Cauchy density is
We can improve this estimate by remembering that -Xi also
has Cauchy distribution. Take Si=-Ti. Remember that Si has
the same distribution as Ti. Then we try (for t>0)
Regression estimates
Suppose we want to compute