Module 8: Probability and Simulation II

We will now return to our estimation example with the uniform distribution:

• Recall that S_n / n is our estimate of the mean.

• For a particular, fixed n (e.g., n=5), S₅ / 4 is a rv with some distribution.

• Consider a slight variation of the sequence: S'_n = (S_n/n - 0.5)

• Clearly, because S_n/n → 0.5), it must be that S'_n → 0

• Let us examine the rv S'_n / n by getting a histogram:

       Algorithm: histogramForS'n_over_n (n)
       Input: n, the fixed number of samples of Xi's.
       1.   for r=1 to numRuns
       2.       total = 0
       3.       for i=1 to n
       4.           x = uniform ()            // The i-th call produces Xi.
       5.           total = total + x         
       6.       endfor
       7.       avg = total / n               // avg = Sn / n
       8.       histogram.add (avg-0.5)       // S'n / n
       9.   endfor
       10.  return histogram
       

• In Java, this looks like: (source file)

          // We'll build a density-histogram in the range [-2,2] with 20 intervals.
          DensityHistogram hist = new DensityHistogram (-2,2,20);
  
          // Number of runs over which to collect samples of S_n/n
          double numRuns = 100000;
  
          for (int r=0; r<numRuns; r++) {
  
              double n = 5;  
              double total = 0;
              for (int i=1; i<=n; i++) {
                  double x = RandTool.uniform ();   // X_i = outcome of i-th call to uniform()
                  total += x;
              }
              double avg = total / n;               // S_n / n
              hist.add (avg-0.5);
          }
  
          hist.display ();
       

  Exercise 22: Download and execute UniformMeanExample3.java). Change the value of n to 50. What do you notice?

• Clearly, because S_n/n - 0.5 → 0, the histograms get narrower and narrower as n increases.
       => Eventually, it becomes a single point.

• On the other hand, consider the sequence S_n by itself.
       => This grows unboundedly.

• Is there something in between, one might ask?

• Consider the sequence S''_n = sqrt(n) * (S_n/n-0.5)

  Exercise 23: Modify UniformMeanExample3.java). to obtain a histogram of S''_n = sqrt(n) * (S_n/n-0.5). Try different values of n.

• What you observe is the Central Limit Theorem in action.

Convergence of rv's to another rv:

• What we saw earlier is that S''_n = sqrt(n) * (S_n/n-0.5) does NOT converge to a number
       => 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₂ + ... + X_n
  • We'll use the symbols μ = E[X₁] and σ = std[X₁] = sqrt(var[X]).
  • Define rv Z_n = (S_n - μ) / σ.
  • Let Z be the standard Normal rv: Z ~ N(0,1)

• The CLT: sqrt(n)*Z_n → Z in distribution.

• Note:
  • The result makes no assumption about the distribution of the X_i's.
  • Z_n = (S_n - μ) / σ is merely a centered (around μ) and scaled (by σ) version of S_n.
  • 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 μ = E[X₁].
  • Accordingly, we obtain an estimate S_n / n, where S_n = (X₁ + ... + 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 μ, we will be interested in whether our estimate "lands" in the interval [μ-δ, μ+δ].
  • Thus, we want to know Pr[estimate lies in [μ-δ, μ+δ]].
         => i.e., Pr[S_n/n ε [μ-δ, μ+δ]].
  • Usually, we want to impose a condition like this:
         => We want Pr[S_n/n ε [μ-δ, μ+δ]] > 0.95.
         => We want to be 95% confident that our estimate is within δ of the true (but unknown) mean μ.
  • For δ, we might want a number like "within 3%"
         => δ = 0.03 μ.

• Now let's pose our two questions in terms of the estimator rv S_n / n:
  • Question 1: How large should n be so that Pr[S_n/n ε [0.97μ, 1.03μ]] > 0.95?
         => 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 δ such that Pr[S_n/n ε [μ-δ, μ+δ]] > 0.95?
         => The interval [μ-δ, μ+δ] 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 S_n/n).

• However, we can approximate the distribution of S_n/n using the CLT and use the Normal distribution to calculate answers to these questions.
       => 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 Pr[S_n/n ε [0.97μ, 1.03μ]] > 0.95?
  • That is, we want that value of n such that
        
  • This is the same as
        
  • Which we can now approximate as
        
  • This is now a calculation to be performed with the Normal distribution.
  • After all, we should be able to compute Pr[a < Z < b] for any interval [a,b].

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

• Here, we will take a further approximation: we will use an estimate of μ and an estimate of σ.

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

• Then, what we need to do is find that value of n such that
      

Estimating μ and σ:

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

• Thus, in the uniform-example, we could use N=1000 samples for the first estimate μ' = S_N/N.

• Now for σ:
  • Recall the definition: var[X] = E[(X-μ)²]
  • 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 μ using μ' = S_N/N.
    2. Step 2: Use the data X₁, ..., X_N and compute the variance estimate
            v = 1/n Σ_i (X_i - μ)²
           
    3. Step 3: Compute
            σ' = sqrt(v).

• Note: in practice it turns out that
       v = 1/(n-1) Σ_i (X_i - μ)²
  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] - &mu')2
       11.   endfor
       12.   σ' = sqrt ( 1/(n-1) varTotal )
       13.   return μ', σ'    
       

• Let's write a program to compute μ' and σ' 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
      

• 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,
         Pr[a < Z < 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,
       2 F_z(a) - 1 > 0.95 => a > 1.96 (From a numerical computation with the CDF).

• Thus, in order that
      
  we need
      
  or
      

Exercise 26: If μ'=0.48 and σ'=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 [μ - fμ, μ + fμ].

• Then, the our earlier result could be stated as:
      

• Thus, the smallest interval is achieved when
      

  Exercise 27: If μ'=0.48 and σ'= 0.27, obtain f when n=500 samples are used.

• Alternatively, we could write the interval as [μ - δ, μ + δ]. => In this case, δ = f μ.

  Exercise 28: Express δ as a function of σ' and n.

Computing the mean, variance and confidence intervals