3
$\begingroup$

I’m trying to better understand several statistical concepts (bootstrapping, central limit theorem, and confidence intervals) by applying them to a binomial distribution (you can think of it as a coin flip, for example).

I’ll explain my expectations/sanity checks, and then hopefully someone can tell me why my expectations don’t match the results I’m getting in my simulation in code.

Procedure:

  1. Given a probability of success (e.g. $p=0.5$), I obtain the results of 1000 bernoulli trials.
  2. I do this 1000 times, which means I get 1000 resamples. From each resample, I compute the sample mean (p=% of positives from the 1000 bernoulli trials), and the sample variance ($n*p*(1-p)$).
  3. I can then construct a “bootstrap distribution” which is the distribution of the sample means.
  4. construct test statistic for $S_n$, standardized and normalized: $$\frac{S_n - p}{\sigma/\sqrt{1000}}$$, where $p$ is the "true" probability of success.
  5. Determine the desired $\alpha$ (type I error rate), and get corresponding $Z$ values: $Z_{\alpha/2}$ and $Z_{1-\alpha/2}$ (for $\alpha=0.05$ I determined these to be $+/- 1.96$.
  6. Determine if the test statistic from step 4) is within the $Z$-values: $Z_{\alpha/2} \le \frac{S_n - p}{\sigma/\sqrt{1000}} \le Z_{1-\alpha/2}$
  7. If so, then we correctly accept the null hypothesis that the probability of success is $p$. If not, then we reject the null hypothesis, and call it a type I error.
  8. Fix $p$ (true probability of success), and $\alpha$ (type I error rate), and repeat steps 1) through 7) 1000 times, and we would expect 50 type I errors.

Expectation (with justification):

  1. I expect the bootstrap distribution to have the following:

    a. mean ($S_n$) = mean of sample means of the resamples (by the law of large numbers).

    b. variance (($\sigma^2)/1000$) = mean of sample variances (by the law of large numbers) divided by 1000 (number of resamples).

    c. Normally distributed (by central limit theorem, see below)

  2. I expect the central limit theorem to apply here because (see “classic CLT” described here): a. The sample means are i.i.d

  3. If the two expectations above are true, then I would expect 95% of the sample means to fall within the confidence interval for the normal distribution parameterized by $S_n$ and $\sigma^2/1000$
  4. I expect the percentage of type I errors to be 5% (50 out of 1000).
    import numpy as np
    from scipy import stats

    def run_test(error_count):
        n=1000 # number of bernoulli trials. one set of bernoulli trials is a resample
        size=1000  # number of resamples
        theta=0.5
        sample_means = []
        sample_variances = []

        # 1) 
        bootstrap_resamples = np.random.binomial(n, theta, size)

        # 2)
        for resample in bootstrap_resamples:
            sample_mean = resample/float(n)
            sample_variance = float(n)*sample_mean*(1-sample_mean)
            sample_means.append(sample_mean)
            sample_variances.append(sample_variance)

        # 3)
        # sample_means (see above) 

        # 4) 
        S_n = np.mean(sample_means)
        sigma_squared = np.mean(sample_variances)
        sigma = np.sqrt(sigma_squared)
        test_statistic  = (S_n - theta)/(sigma/np.sqrt(1000))

        # 5) 
        alpha=0.05
        if not stats.norm.ppf(alpha/2.0) <  test_statistic < stats.norm.ppf(1-(alpha/2.0)):
            error_count +=1
        return error_count


    error_count = 0
    for i in range(1000):
        error_count = run_test(error_count)

    print error_count

The above code consistently returns 0 errors. In other words, print error_count returns 0, when I would expect it to return approximately 50.

In most cases it might be good to have 0 errors, but I'm trying to validate that the $\alpha$ level influences the number of errors, and I'm not seeing that here.

I have a suspicion that the test statistic is computed incorrectly, but I can't figure out what it is...Do you perhaps see an issue in my computation?

Thank you.

$\endgroup$
16
  • 1
    $\begingroup$ 1. "I always thought that the central limit theorem meant that the sample means would be distributed by a standard normal." -- this is wrong in several ways at once. It might be worth reading about various forms of the CLT on Wikipedia (but perhaps focus on the classical CLT for simplicity). 2. Posting code without first clearly saying what the code is implementing is problematic because it conflates errors of understanding with errors of implementation. It's much better to resolve your errors of understanding first (by explaining what you're intending to implement), without code. $\endgroup$
    – Glen_b
    Feb 20, 2020 at 5:30
  • 1
    $\begingroup$ @Glen_b-ReinstateMonica, I did what you suggested, and I looked at classical CLT from en.wikipedia.org/wiki/Central_limit_theorem. I hope the question looks much nicer now. Thanks. $\endgroup$
    – makansij
    Feb 21, 2020 at 21:05
  • 3
    $\begingroup$ This looks fine to me--except you do not describe a bootstrap. It's close to a parametric bootstrap, but subtly different. The parametric bootstrap would be executed by first generating a sample of 1000 iid Bernoulli(0.5) variates and setting $p$ to equal their mean. Then you would resume at step 1 using this $p$ rather than a value of $0.5.$ The point is that in the bootstrap you never know the true mean: you are using the sample estimate as a surrogate. $\endgroup$
    – whuber
    Feb 21, 2020 at 21:39
  • 1
    $\begingroup$ My initial hesitation was simply trying to guess what you were really doing (rather than trying to figure out what your code did). $\endgroup$
    – Glen_b
    Feb 23, 2020 at 1:09
  • 2
    $\begingroup$ Could you explain what "giving me 0 errors" means? One might interpret that as "giving me zero errors; i.e.,, no errors," which ideally is what you want, right? $\endgroup$
    – whuber
    Feb 23, 2020 at 21:43

1 Answer 1

2
$\begingroup$

Let's first describe with mathematical notation what happens in the simulation (when dividing by n) and then lay out the statistics that would lead one to expect the result.

bootstrap_resamples = np.random.binomial(n, theta, size)

Generates $S=1000$ binomial random variables, $Y_s$ each consisting of the sum of $N=1000$ Bernoulli trials with success probability $\theta$, $Y_s = \sum_j^{N} X_{s, j}$

sample_mean = resample/float(n)
sample_variance = sample_mean*(1-sample_mean) / n

Calculates $\hat{\theta}_s = \frac{1}{N} \sum_j^{N} X_{s, j}$ and $\hat{\text{v}}_s = \frac{\hat{\theta_s}(1-\hat{\theta_s})}{N}$

  S_n = np.mean(sample_means)
  sigma_squared = np.mean(sample_variances)
  sigma = np.sqrt(sigma_squared)
  test_statistic  = (S_n - theta)/(sigma/np.sqrt(1000))

Calculates $\hat{\theta} = \frac{1}{S} \sum_s^S \hat{\theta}_s$ and $\hat{\text{v}} = \frac{1}{S} \sum_s^S \hat{\text{v}_s}$ then $\text{stat} = \frac{ \hat{\theta} - \theta}{ \sqrt{ S^{-1} \hat{\text{v}} } }$

       if not stats.norm.ppf(alpha/2.0) <  test_statistic < stats.norm.ppf(1-(alpha/2.0)):
        error_count +=1

Finally, compute the indicator variable $Q = \text{stat} \not \in (-1.96, 1.96)$

Now, what is the sampling distribution of $Q$? If $\text{stat}$ is approximately standard normal, then it approximately equals $5\%$. This is in fact the case because a central limit theorem applies.

Looking at $\hat{\theta} = \frac{1}{S \cdot N} \sum_s \sum_j X_{s,j}$ that is just an average of $S\cdot N$ Bernoulli variables. If we subtract its mean divide by its standard deviation, it will have mean zero and variance 0. It consists of a sum of iid random variables, so will be approximately normal in large samples and be within (-1.96, 1.96) $95\%$ of the time. We have $E[ \hat{\theta} ] = \theta$ and $Var( \hat{\theta} ) = \text{V} = \frac{\theta(1-\theta)}{S\cdot N} $.

So $$\frac{ \hat{\theta} - \theta }{ \sqrt{\text{V}}}$$ is approximately standard normal by a CLT. In the simulation, the approximation will be really good because the sum is over a million random variables. If you want to focus only on understanding the CLT, divide by $v$ in the simulation. And you could have just looked at the mean of 1000 bernoulli trials instead of the mean over 1000 means.

However, the denominator of $\text{stat}$ isn't $\sqrt{V}$ and yet the simulation still yields $Q$ at around $5\%$ anyways. Why? Because the denominator of $\text{stat}$ converges in probability to the root of $V$ by a law of large numbers. And there's another theorem (https://en.wikipedia.org/wiki/Slutsky%27s_theorem) that says asymptotically, that's as good as if it it actually was $V$. To see the denominator works out: The denominator of $\text{stat}$ is the root of $\frac{1}{S \cdot S} \sum_s^S \hat{\text{v}_s} = \frac{1}{S \cdot S} \sum_s^S \frac{\hat{\theta_s}(1-\hat{\theta_s})}{N} = \frac{1}{S \cdot N} \left( \frac{1}{S} \sum_s \hat{\theta_s} - \frac{1}{S} \sum_s \hat{\theta_s}^2 \right)$. The quantities inside the brackets will converge to $\theta$ and $\theta^2$ by a law of large numbers, so the denominator will converge to the root of $V$ and Slutsky's theorem kicks in.

$\endgroup$
3
  • $\begingroup$ This is extremely helpful. How did you know to divide by $N$ in the formula for estimating the variance of $S$? I could understand if that were the formula for the variance of a bernoulli distribution. And I understand that the mean of a Bernoulli distribution is obviously just the probability of success, which in our case is taken to be $\hat \theta_s$, as you have indicated. But the corresponding variance $v_S$ of a Bernoulli distribution is $v_S = \hat \theta_s ( 1 - \hat \theta_s )$, but you have $v_S = \frac{\hat \theta_s ( 1 - \hat \theta_s )}{N}$. How did you know to divide by $N$? $\endgroup$
    – makansij
    Mar 1, 2020 at 4:39
  • 1
    $\begingroup$ @Hunle: The variance is only needed to scale the statistic so that its variance is $1$ (otherwise, its distribution cannot be compared to a standard normal random variable which has variance 1). The statistic here is a sample mean $\frac{1}{N} \sum_i X_i$, and the variance of a sample mean of uncorrelated random variables is s $\frac{\sigma^2}{N}$ $\endgroup$
    – CloseToC
    Mar 3, 2020 at 16:19
  • $\begingroup$ I see, so in this case, we’re computing the variance of a bernoulli distribution, and then scaling it to fit the standard normal as desired. The independent random variables are the number of successes from each sample of bernoulli trials (as you denote it, $X_{s,j}$. The sample mean is obviously $\hat \theta_S$, as you indicated, and since these are means of bernoulli trials, the variance is $\hat \theta_S(1-\hat \theta_S)$. Then, we shift and scale the distribution of the sample mean by first subtracting the true mean 0.5, and scaling by square root of sample variance, divided by $N$. $\endgroup$
    – makansij
    Mar 3, 2020 at 21:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.