Module 9: Probability and Simulation II

Convergence of rv's to another rv:

• What we saw earlier is that \(S^{''}_n = \sqrt{n} (\frac{S_n}{n}-0.5)\) does NOT converge to a number
       \(\Rightarrow\) It converges to another rv (whose pdf happens to be bell-shaped).

• Thus, we now have another type of limit: convergence in distribution.

• Formal definition: the sequence \(Y_n\) converges to rv \(Y\) if the CDFs of \(Y_n\) converge to the CDF of \(Y\).

• Note: this involves the idea of limits of functions (the CDFs).

• Thus, the Central Limit Theorem is really about converging to another rv (or, equivalently, about converging to another distribution).

The Central Limit Theorem:

• First, the background:
  • Let \(X_i\) be the outcome of the \(i\)-th independent trial of an experiment.
  • Let \(S_n = X_1 + X_2 \ldots + X_n\)
  • We'll use the symbols \(\mu = \expt{X_1}\) and \(\sigma = \mbox{std}[X_1] = \sqrt{\mbox{var}[X]}\).
  • Define rv $$ Z_n \eql \frac{(\frac{S_n}{n} - \mu)}{\sigma} $$
  • Let \(Z\) be the standard Normal rv: \(Z \simm{N}(0,1)\)

• The CLT: \(\sqrt{n}Z_n \to Z\) in distribution.

• Note:
  • The result makes no assumption about the distribution of the \(X_i\)'s.
  • The variable $$ Z_n \eql \frac{(\frac{S_n}{n} - \mu)}{\sigma} $$ is merely a centered version of \(S_n\), centered around \(\mu\) and scaled by \(\sigma\).
  • If we didn't scale, it would still work, but it would converge to a non-standard Normal.

Using the CLT to determine estimation quality:

• Recall our estimation problem:
  • We are repeating an experiment and each repetition provides a sample \(X_i\).
  • We are interested in estimating the true value of \(\mu = \expt{X_1}\).
  • Accordingly, we obtain an estimate \(\frac{S_n}{n}\), where \(S_n = (X_1 + \ldots + X_n)\).

• There are two key questions related to estimation quality:
  • If we are at liberty to choose \(n\), how large should \(n\) be?
  • If not, then, for fixed \(n\), how accurate is the estimate likely to be?

• We need some notion of accuracy here:
  • Since we are estimating \(\mu\), we will be interested in whether our estimate "lands" in the interval \([\mu-\delta, \mu+\delta]\).
  • Thus, we want to know \(\prob{\mbox{estimate lies in } [\mu-\delta, \mu+\delta]}\), i.e. $$ \prob{\frac{S_n}{n} \in [\mu-\delta, \mu+\delta]} $$
  • Usually, we want to impose a condition like this:
         \(\Rightarrow\) We want $$ \prob{\frac{S_n}{n} \in [\mu-\delta, \mu+\delta]} \gt 0.95 $$
         \(\Rightarrow\) We want to be 95% confident that our estimate is within \(\delta\) of the true (but unknown) mean \(\mu\).
  • For \(\delta\), we might want a number like "within 3%"
         \(\Rightarrow\) \(\delta = 0.03 \mu\).

• Now let's pose our two questions in terms of the estimator rv \(\frac{S_n}{n}\):
  • Question 1: How large should \(n\) be so that $$ \prob{\frac{S_n}{n} \in [0.97\mu, 1.03\mu]} \gt 0.95? $$
         \(\Rightarrow\) i.e., What is \(n\) so that we're 95% sure that our estimate will be within 3% of the mean?
  • Question 2: For fixed \(n\), what is the value of \(\delta\) such that $$ \prob{\frac{S_n}{n} \in [\mu-\delta, \mu+\delta]} \gt 0.95? $$
         \(\Rightarrow\) The interval \([\mu-\delta, \mu+\delta]\) in this case is called a confidence interval.

• We can't really compute the answer exactly because we don't know the distribution of \(X_i\)'s (much less the distribution of \(\frac{S_n}{n}\)).

• However, we can approximate the distribution of \(\frac{S_n}{n}\) using the CLT and use the Normal distribution to calculate answers to these questions.
       \(\Rightarrow\) This calculation is not straightforward, so let's take a closer look.

Question 1: how many samples are needed?

• Suppose we want to know how large \(n\) should be so that $$ \prob{\frac{S_n}{n} \in [0.97\mu, 1.03\mu]} \gt 0.95? $$
  • That is, we want that value of \(n\) such that $$ \prob{\mu-0.03\mu \; \lt \; \frac{S_n}{n} \; \lt \; \mu+0.03\mu} \;\; \gt \;\; 0.95 $$
  • This is equivalent to (subtracting off \(\mu\)): $$ \prob{-0.03\mu \; \lt \; \frac{S_n}{n} -\mu \; \lt \; 0.03\mu} \;\; \gt \;\; 0.95 $$
  • Which is the same as (dividing by \(\sigma\)): $$ \prob{- \frac{0.03\mu}{\sigma} \; \lt \; \frac{\frac{S_n}{n} - \mu}{\sigma} \; \lt \; \frac{0.03\mu}{\sigma} } \;\; \gt \;\; 0.95 $$
  • Which is the same as (multiplying by \(\sqrt{n}\)): $$ \prob{- \frac{\sqrt{n}0.03\mu}{\sigma} \; \lt \; \sqrt{n}\frac{\frac{S_n}{n} - \mu}{\sigma} \; \lt \; \frac{\sqrt{n}0.03\mu}{\sigma} } \;\; \gt \;\; 0.95 $$
  • Which we can now approximate as $$ \prob{- \frac{\sqrt{n}0.03\mu}{\sigma} \; \lt \; Z \; \lt \; \frac{\sqrt{n}0.03\mu}{\sigma} } \;\; \gt \;\; 0.95 $$ since $$ \sqrt{n}\frac{\frac{S_n}{n} - \mu}{\sigma} \; \approx \; Z $$
  • This is now a calculation to be performed with the Normal distribution.

  • After all, we are able to compute \(\prob{a \lt Z \lt b}\) for any interval \([a,b]\) since we know the distribution of \(Z\).

• However, there is a small problem: doesn't this imply that we have to already know \(\mu\) and \(\sigma\)?
  • Recall, we are trying to estimate the number \(\mu\)
         \(\Rightarrow\) Yet, the above calculation seems to need it.
  • Similarly, not knowing anything about \(X_i\), how can we compute its standard deviation \(\sigma\)?

• Here, we will take a further approximation: we will use an estimate of \(\mu\) and an estimate of \(\sigma\).

• Let's say we somehow (we'll talk about this shortly) have estimates:
  • Let \(\mu'\) denote our estimate of \(\mu\).
  • Let \(\sigma'\) denote our estimate of \(\sigma\).

• Then, what we need to do is find that value of \(n\) such that $$ \prob{- \frac{\sqrt{n}0.03\mu'}{\sigma'} \; \lt \; Z \; \lt \; \frac{\sqrt{n}0.03\mu'}{\sigma'} } \;\; \gt \;\; 0.95 $$

Estimating \(\mu\) and \(\sigma\):

• First, let's tackle \(\mu\):
  • We know how to estimate \(\mu\): by obtaining \(\frac{S_n}{n}\).
  • There seems to be some circular logic here: to know how many samples we need, we first have to estimate \(\mu\).
  • The way out is to use as large a number as you can initially: pick some large \(N\) (e.g., 1000) and let
         \(\mu^\prime = \frac{S_n}{N}n\)

• Thus, in the uniform-example, we could use \(N=1000\) samples for the first estimate \(\mu^\prime = \frac{S_n}{N}\).

• Now for \(\sigma\):
  • Recall the definition: \(\mbox{var}[X] = \expt{(X-\mu)^2}\)
  • To estimate this, we need to compute from data the average value of this quantity.
  • We'll do this in two steps:
    1. Step 1: estimate \(\mu\) using \(\mu^\prime = \frac{S_n}{N}\).
    2. Step 2: Use the data \(X_1, X_2, \ldots, X_n\) and compute the variance estimate $$ v \eql \frac{1}{n} \sum_i (X_i - \mu^\prime)^2 $$
    3. Step 3: Compute
            \(\sigma^\prime = \sqrt{v}\).

• Note: in practice it turns out that $$ v \eql \frac{1}{n-1} \sum_i (X_i - \mu^\prime)^2 $$ is a better estimate of the variance.

• Thus, in an application

       Algorithm: estimateMuSigma (n)
       Input: n, samples to use in this early estimate
       1.    muTotal = 0
       2.    for i=1 to n
       3.        x = outcome from experiment
       4.        muTotal = muTotal + x
       5.        data[i] = x                           // Store the data for now.
       6.    endfor
       7.    μ' = muTotal / n
             // Next compute the standard deviation.
       8.    varTotal = 0
       9.    for i=1 to n
       10.       varTotal = varTotal + (data[i] - μ')^2
       11.   endfor
       12.   σ' = sqrt ( 1/(n-1) varTotal )
       13.   return μ', σ'    
       

• Let's write a program to compute \(\mu^\prime\) and \(\sigma^\prime\) for our uniform-example: (source file)

          double n = 100;                
  
          // First estimate the mean.
          double total = 0;
          ArrayList<Double> data = new ArrayList<Double> ();
          for (int i=0; i<n; i++) {
              double x = RandTool.uniform ();   
              total += x;
              data.add (x);                       // Store for later use in variance estimate.
          }
          double mu_p = total / n;       
          System.out.println ("Mean estimate: " + mu_p);   
          
          // Now the variance.
          double varTotal = 0;
          for (double x: data) {
              varTotal += (x - mu_p) * (x - mu_p);
          }
          double sigma_p = Math.sqrt ( (1.0/(n-1) * varTotal) );
          System.out.println ("Std-dev estimate: " + sigma_p);
  
          // One can now use mu_p and sigma_p in Normal-distribution formulas.
        

Exercise 24: Execute UniformVarianceExample.java and obtain the two estimates.

Finally, let's go back to our "how many samples" question:

• Recall that we want that value of \(n\) such that $$ \prob{- \frac{\sqrt{n}0.03\mu'}{\sigma'} \; \lt \; Z \; \lt \; \frac{\sqrt{n}0.03\mu'}{\sigma'} } \;\; \gt \;\; 0.95 $$

• At this point, we could simply look this up in a book or write code to estimate this quantity.

• However, let's take this forward a few steps:
  • Let \(F_z\) denote the CDF of the standard Normal rv.
  • Then, for any value \(a\),
         \(\prob{a \lt Z \lt a} = 2 F_z(a) - 1\)

  Exercise 25: This last result is actually not difficult to prove. Draw a picture of the Normal pdf and see if you can prove this.

• Then, one can use the CDF (numerically) to show that
       \(2 F_z(a) - 1 \gt 0.95\)
  when
       \(a \gt 1.96\) (From a numerical computation with the CDF).

• Thus, in order that $$ \prob{- \frac{\sqrt{n}0.03\mu'}{\sigma'} \; \lt \; Z \; \lt \; \frac{\sqrt{n}0.03\mu'}{\sigma'} } \;\; \gt \;\; 0.95 $$
  we need $$ \frac{\sqrt{n}0.03\mu'}{\sigma'} \;\; \gt \;\; 1.96 $$
  or $$ n \; \gt \; \left( \frac{1.96\sigma'}{0.03\mu'} \right)^2 $$

Exercise 26: If \(\mu'=0.48\) and \(\sigma'=0.27\) and obtain the number of samples needed.

Lastly, the second question: for a fixed \(n\) what is the width of the confidence interval?

• Suppose \(f\) is the fraction so that the interval is \([\mu - f\mu, \mu + f\mu]\).

• Then, the our earlier result could be stated as: $$ f \; \gt \; \frac{1.96 \sigma'}{\sqrt{n} \mu'} $$

• Thus, the smallest interval is achieved when $$ f \eql \frac{1.96 \sigma'}{\sqrt{n} \mu'} $$

  Exercise 27: If \(\mu'=0.48\) and \(\sigma'= 0.27\), obtain \(f\) when \(n=500\) samples are used.

• Alternatively, we could write the interval as \([\mu - \delta, \mu + \delta]\). \(\Rightarrow\) In this case, \(\delta = f \mu\).

  Exercise 28: Express \(\delta\) as a function of \(\sigma'\) and \(n\).

9.10   Computing the mean, variance and confidence intervals

Let's write a simple statistics object to track basic statistics in a simulation:

• First, note that the sample mean is easy to track:

  public class Stats {
  
      double numSamples;           // Count # samples received so far.
      double sampleSum;            // Accumulated sum.
      double sampleMean;           // Mean.
  
      public void add (double x)
      {
          numSamples ++;
          sampleSum += x;
          sampleMean = sampleSum / numSamples;
      }
  
      public double getMean ()
      {
          return sampleMean;
      }
      
  }
  
       

• However, recall our definition of variance:
  • We first estimate the mean \(\mu' = \frac{S_n}{N}\).
  • Then, estimate the variance $$ v \eql \frac{1}{n} \sum_i (X_i - \mu^\prime)^2 $$
  Doesn't this mean that we need to store all the samples?

• It turns out that one can express the variance in a recurrence relation as follows:
  • Let \(\mu_n\) = estimate of the mean using the first \(n\) samples.
  • Let \(v_n\) = estimate of the variance using the first \(n\) samples.
  Then, one can show that $$ v_n \eql \left( \frac{n-2}{n-1} \right) v_{n-1} \; + \; n \left( \mu_n - \mu_{n-1} \right)^2 $$ Which means there's no need to store the data values.

• Accordingly, to track the variance:

      public void add (double x)
      {
          numSamples ++;
          sampleSum += x;
          double oldMean = sampleMean;
          sampleMean = sampleSum / numSamples;
          if (numSamples <= 1) {
              // This should be treated as an improper value.
              sampleVar = 0;
          }
          else {
              // Note: recurrent formula for variance \(\Rightarrow\)  no need to store samples.
              sampleVar = (1 - 1.0/(double)(numSamples-1)) * sampleVar 
                           + numSamples * (sampleMean - oldMean)*(sampleMean - oldMean);
          }
      }
       

• Finally, let's add code to compute confidence intervals:

          // We're using 95% confidence here. 
          delta = 1.96 * Math.sqrt (sampleVar/numSamples);