The mathematics of the long run

(Recurring ideas)

Andrew Pua

2025-10-20

Example of random processes

  1. Suppose that on one night at a certain hospital, four mothers give birth to baby boys.
  2. As a very sick joke, the hospital staff decides to return the babies to their mothers completely at random.
  3. How do you use R to implement the sick joke and check how many babies were returned to the correct mother?
  4. Repeat the R code many times. What do you notice?

Exploration of a random process

  1. What if you repeat this for 100 times? Is there a more systematic way to do this?

  2. Can you summarize the information?

    • What seems to be the most likely value for the number of matches? Least likely?
    • What is the average number of matches?
  3. How many repetitions do you need to get some form of stability?

  4. Do you notice something from your summary?

Probability

  1. The probability of an outcome refers to how often the outcome would occur in the long run if a random process were repeated over and over under identical conditions.

  2. Probability is a long run concept. The empirical counterpart is the relative frequency.

  3. What is your empirical estimate of the probability of all four mothers getting the correct baby?

  4. Can you create R code which will enable you to determine how the average number of matches will behave as you increase the number of repetitions?

Language

  1. Sample Space: List all possible outcomes systematically.
  2. Use given conditions to assign numbers between 0 and 1 to each outcome.
  3. Probability Function: But these numbers cannot be negative and should sum up to 1 over all outcomes.
  4. Random assignment means all possible outcomes in the sample space are equally likely.
  1. What is the probability that all four mothers get the correct baby?

  2. The set of outcomes which correspond to “all four mothers get the correct baby” is called an event.

  3. Event: A collection of outcomes from the sample space

  4. For each outcome in the sample space, determine how many mothers received the right baby.

  1. You created a random variable. It is usually denoted by a capital letter, say \(X\).

    • A random variable takes each outcome in the sample space and assigns it a numerical value.
    • A random variable is discrete if you can list out each individual value that can be observed for the random variable.
  2. After determining how many mothers received the right baby, determine the probability that none of the mothers receives the correct baby.

  3. Repeat this for all possible value for the number of correct matches.

  1. You have created the probability distribution of this discrete random variable, the list of all possible values and the probability assigned to each value.
  2. Sketch a graph of this probability distribution.
  3. How do these probabilities compare to the relative frequencies earlier?
  4. Determine the probability that at least one mother gets the correct baby. How does this relate to the probability that no mothers will get the right baby?
  1. Recall that the average number of matches will tend to stabilize too.
  2. The expected value of a random variable is the long run average of the random variable.
  3. For a discrete random variable, it is calculated by taking a probability-weighted average. Let \(x_1,\ldots,x_K\) denote the list of possible outcomes of a random variable \(X\). Then, we have \[\mathbb{E}\left(X\right)=\sum_{i=1}^K x_i\mathbb{P}\left(X=x_i\right)\]
  1. How does this expected value compare to the long run average value from the Monte Carlo simulation?

  2. Does the expected value necessarily equal the most likely value? Does it even have to be a possible value?

  3. Another property of a random variable is its variance. This measures how variable or how spread out the outcomes of the random variable will be relative to the expected value.

  1. For a discrete random variable, it also has the form of a probability-weighted average. Let \(x_1,\ldots,x_K\) denote the list of possible outcomes of a random variable \(X\). Let \(\mathbb{E}\left(X\right)\) be the expected value. Then, we have \[\mathsf{Var}\left(X\right)=\sum_{i=1}^K \left(x_i-\mathbb{E}\left(X\right)\right)^2\mathbb{P}\left(X=x_i\right)\]

  2. The square root of the variance of a random variable is the standard deviation of a random variable, i.e. \(\mathsf{SD}\left(X\right)=\sqrt{\mathsf{Var}\left(X\right)}\).

  1. What is the difference between the standard deviation you calculate using data and the standard deviation of a random variable?

  2. How will all calculations change if there are 8 babies? \(B\) babies?

Skills, mysteries, and connections

  1. Need to know counting methods, instead of listing all possible outcomes
  2. Need to know how to do a Monte Carlo simulation in case you don’t know how to count
  3. Need to be aware of other assignment mechanisms for probabilities
  4. What happens if we are not in the discrete case?
  5. How is everything we have seen so far connected to what we did for hypothesis testing, for example?

Revisiting Harley

  1. A binomial random variable counts the number of “successes” in a random process with the following properties:

    • Each trial results in “success” or “failure”.
    • The trials are independent: The outcome of one trial does not change the probability of success on the next trial.
    • The probability of success, \(\pi\), is constant across the trials.
    • There are a fixed number of trials, \(n\).
  1. If a random variable \(X\) counts the number of successes, then we can say \(X\sim \mathsf{Bin}\left(n, \pi\right)\), or \(X\) has a binomial distribution with parameters \(n\) and \(\pi\).

  2. It can be shown that a binomial random variable has the following properties:

    • \(\mathbb{P}\left(X=x\right)=\dfrac{n!}{x!\left(n-x\right)!} \pi^x (1-\pi)^{n-x},\ x=0,1,\ldots,n\)
    • \(\mathbb{E}\left(X\right)=n\pi\)
    • \(\mathsf{Var}\left(X\right)=n\pi\left(1-\pi\right)\)
  1. The probability distribution of \(X\) can be obtained by applying two probability rules:

    • When the trials are independent (the outcome on one trial does not affect the probability of success on another trial), the probability of several outcomes occurring simultaneously is the product of the probabilities of the individual outcomes.
    • The probability of one of several mutually exclusive events (events that can’t occur simultaneously) occurring is the sum of the probabilities of the individual outcomes.
  1. With Harley, we were interested in how likely is it to see a relative frequency of 0.9 (9 out of 10 times Harley chose the correct cup) if Harley did not understand human gestures.

    • Define \(X\) to be the number of times Harley chose the correct cup.
    • Under the conditions of the experiment, \(X\sim \mathsf{Bin}\left(n,\pi\right)\) seems plausible. Here \(n=10\).
    • But what is \(\pi\)?
  1. Under the “ordinary” situation that Harley does not understand human gestures, \(\pi\leq 0.5\).

  2. Under the “out-of-ordinary” situation that Harley does understand human gestures, \(\pi > 0.5\).

  3. We want the probability that \(X\geq 9\) under the “ordinary” situation. This means we want \(\mathbb{P}\left(X\geq 9\right)\) under \(X\sim \mathsf{Bin}\left(10,0.5\right)\).

  4. Find this value and compare with our simulation before.

Revisiting left-handed scientists

  1. About 13% of the population is left-handed. A biologist suspects that the scientific community is not like the general population in terms of handedness. He asks 200 scientists and finds that 36, are left-handed. Does this data support the biologist’s theory?

  2. Define a random variable \(X\) be the number of scientists who are left-handed.

  3. How plausible do you think is the assumption that \(X\) is a binomial random variable?

  1. If we are convinced of its plausibility, then we have \(X\sim \mathsf{Bin}\left(200, 0.13\right)\).

  2. But this was a situation where the “out-of-ordinary” claim could go in two directions.

  3. It was relatively straightforward to compute \(\mathbb{P}\left(X\geq 36\right)\) when \(X\sim\mathsf{Bin}\left(200,0.13\right)\). Find the value.

  4. But the earlier calculation was only for the right tail probability How about the left tail probability?

  5. Last time, I mentioned that we could double the \(p\)-value obtained. But this works only under certain conditions.

  1. One approach is to exploit the fact that \(X\) is a binomial random variable.

    • First, find \(\mathbb{E}\left(X\right)\) under the “ordinary” situation.
    • Second, find the \(x\) such that \(\mathbb{P}\left(X\leq x\right) \leq \mathbb{P}\left(X\geq 36\right)\). This is done by trial and error and taking the \(\mathbb{E}\left(X\right)\) under the “ordinary” situation.
    • The two-tailed \(p\)-value is then \(\mathbb{P}\left(X\leq x\right)+\mathbb{P}\left(X\geq 36\right)\).
  2. But this can be difficult to calculate by hand. So, use R.

# expected value of X under the "ordinary" situation
ev <- 200*0.13
# right tail probability
rtp <- pbinom(35, 200, 0.13, lower.tail = FALSE)
rtp
[1] 0.02667494
# possible values of P(X<=x) for x <= ev 
probs <- pbinom(0:ev, 200, 0.13, lower.tail = TRUE)
probs
 [1] 8.014022e-13 2.475135e-11 3.808348e-10 3.892554e-09 2.973599e-08
 [6] 1.811131e-07 9.162493e-07 3.960606e-06 1.493516e-05 4.991919e-05
[11] 1.497644e-04 4.074632e-04 1.013944e-03 2.324500e-03 4.940232e-03
[16] 9.786852e-03 1.816050e-02 3.170328e-02 5.227690e-02 8.172467e-02
[21] 1.215469e-01 1.725506e-01 2.345600e-01 3.062689e-01 3.852929e-01
[26] 4.684226e-01 5.520300e-01
# which x is below the threshold
max(which(probs <= rtp))
[1] 17
# left-tail probability
probs[max(which(probs <= rtp))]
[1] 0.0181605

Effect of changing \(n\) and \(\pi\)

  1. You will see pictures of the probability distribution of \(X\sim\mathsf{Bin}\left(n,\pi\right)\).

  2. \(n\) could be 10, 100, 1000. \(\pi\) could be 0.5 or 0.05.

  3. You could get a sense of how the shape of the probability distribution is affected by \(n\) and \(\pi\).

  4. The pictures give you a sense of when doubling the \(p\)-value works.

\(n\) increases, \(\pi\) held constant

pi <- 0.5; n <- 10; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")
pi <- 0.5; n <- 100; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")
pi <- 0.5; n <- 1000; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")

\(n\) increases, \(\pi\) held constant

pi <- 0.05; n <- 10; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")
pi <- 0.05; n <- 100; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")
pi <- 0.05; n <- 1000; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")

\(n\) fixed, \(\pi\) changing

pi <- 0.05; n <- 100; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")
pi <- 0.1; n <- 100; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")
pi <- 0.15; n <- 100; x <- seq(0, n, 1)
plot(x, dbinom(x, n, pi), type = "h")

Estimating \(\pi\) when \(n\) is given

  1. You already have an intuitive feel that we can estimate \(\pi\) as the relative frequency of success in a sample.

  2. From a procedural standpoint, we are dividing a binomial random variable \(X\) by \(n\), i.e. a point estimator of \(\pi\) is \(\widehat{\pi}=X/n\).

  3. In hypothesis testing, we have a claim about \(\pi\) and we have a point estimate – the observed value of \(\widehat{\pi}\).

  4. Another common situation is that we do not know anything about \(\pi\).

  1. You might be asked to provide a set of plausible values of the \(\pi\) based on the the observed value of \(\widehat{\pi}\).

  2. This leads to the notion of a confidence interval.

  3. We can employ a “trial-and-error” approach to determine which values of \(\pi\) appear plausible based on what we observed in the sample.

  4. This involves testing different values of \(\pi\) and seeing whether the corresponding two-sided \(p\)-value is larger than some pre-specified cut-off \(\alpha\).

  1. You have to set the value of \(\alpha\). In practice \(\alpha=0.05\) is used frequently.

  2. Consider our left-handed scientists example again. We want a set of plausible values of \(\pi\).

    • Adjust the value of \(\pi\) so that \(\mathbb{P}\left(\widehat{\pi}\geq 0.18\right) \approx 0.025\).
    • Adjust the value of \(\pi\) so that \(\mathbb{P}\left(\widehat{\pi}\leq 0.18\right) \approx 0.025\).
  3. Another way to say this is as follows:

    • Adjust the value of \(\pi\) so that \(\mathbb{P}\left(X\geq 36\right) \approx 0.025\).
    • Adjust the value of \(\pi\) so that \(\mathbb{P}\left(X\leq 36\right) \approx 0.025\).
# Start trying 0.11 to 0.15 as plausible values of pi given what we see in the sample
trial.values <- seq(0.11, 0.15, 0.01)
# Get a sense of how to narrow the location
pbinom(35, 200, trial.values, lower.tail = FALSE)
[1] 0.002125608 0.008566697 0.026674944 0.066761037 0.138732538
# Try a narrower range
trial.values <- seq(0.12, 0.13, 0.0005)
pbinom(35, 200, trial.values, lower.tail = FALSE)
 [1] 0.008566697 0.009119709 0.009702233 0.010315443 0.010960538 0.011638739
 [7] 0.012351290 0.013099456 0.013884523 0.014707797 0.015570602 0.016474279
[13] 0.017420186 0.018409696 0.019444196 0.020525083 0.021653767 0.022831666
[19] 0.024060207 0.025340821 0.026674944
# More systematic way
f1 <- function(x) {
  pbinom(35, 200, x, lower.tail = FALSE) - 0.025
}
uniroot(f1, c(0.12,0.13))
$root
[1] 0.1293614

$f.root
[1] -1.95471e-05

$iter
[1] 3

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05
# More systematic way, the other tail
f2 <- function(x) {
  pbinom(36, 200, x, lower.tail = TRUE) - 0.025
}
uniroot(f2, c(0.19,0.50))
$root
[1] 0.2403667

$f.root
[1] 2.700976e-05

$iter
[1] 9

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05
  1. Therefore, a 95% confidence interval for \(\pi\) given the sample result is
uniroot(f1, c(0.05,0.15))$root
[1] 0.129368
uniroot(f2, c(0.19,0.50))$root
[1] 0.2403667
  1. The approach here is not very convenient, but the approach here illustrates this idea directly.

  2. There are formulas developed so that you can avoid what I have shown you.

  3. Next, you see a visualization.

c1 <- uniroot(f1, c(0.05,0.15))$root 
c2 <- uniroot(f2, c(0.19,0.50))$root
n <- 200; x <- seq(0, n/3, 1) # to make things easier to see
par(mfrow=c(2,1))
plot(x, dbinom(x, n, c1), type = "h")
abline(v = 36, col = "red")
plot(x, dbinom(x, n, c2), type = "h")
abline(v = 36, col = "red")
  1. The final situation is when you are content with giving a point estimate of \(\pi\) and express uncertainty about the fact that one’s point estimate is subject to sampling variation.

  2. You already know that \(\widehat{\pi}\) or the relative frequency of “success” is a good way to estimate \(\pi\), at least when \(n\) is large enough.

  3. To express the fact that \(\widehat{\pi}\) is subject to variation, we need not look further than reporting the SD of \(\widehat{\pi}\).

  1. We now show that \[\mathsf{SD}\left(\widehat{\pi}\right) = \sqrt{\frac{\pi\left(1-\pi\right)}{n}}\]

  2. We need a couple of properties of expected values and variances. Let \(c\) be a constant and \(X\) be any random variable. You can prove that:

    • \(\mathbb{E}\left(cX\right)=c\mathbb{E}\left(X\right)\)
    • \(\mathsf{Var}\left(cX\right)=c^2\mathsf{Var}\left(X\right)\)
  1. Taking \(X\sim \mathsf{Bin}\left(n,\pi\right)\) and noting that \(\widehat{\pi}=X/n\), we let \(c=1/n\) to obtain

    • \(\mathbb{E}\left(\widehat{\pi}\right)= \dfrac{1}{n}\times n\pi = \pi\)
    • \(\mathsf{Var}\left(\widehat{\pi}\right)=\dfrac{1}{n^2}\times n\pi\left(1-\pi\right)=\dfrac{\pi\left(1-\pi\right)}{n}\)
    • Take square roots and you obtain the result.
  1. \(\mathsf{SD}\left(\widehat{\pi}\right)\) is sometimes called the standard error of \(\widehat{\pi}\).

    • Why introduce a new term?
    • Standard error applies to hypothetical proportions which could have been observed.
  2. The problem is that the standard error SE depends on an unknown \(\pi\). So we have to plug-in a “good” estimate for \(\pi\).

  3. Therefore, the estimated standard error \(\widehat{\mathsf{SE}}\) is given by \[\widehat{\mathsf{SE}}=\sqrt{\frac{\widehat{\pi}\left(1-\widehat{\pi}\right)}{n}}\]

Seeing the results through a simulation

  1. You will now see a Monte Carlo simulation of the behavior of \(\widehat{\pi}\).

    • Must generate simulated data via rbinom().
    • Must think about the underlying parameters to generate simulated data.
  2. Settings:

    • \(n\) could be 8, 32, 128, …
    • \(\pi=0.5\)
dist.mean <- function(n, pi)
{
  a <- rbinom(10^4, n, pi)
  sample.means <- a/n
  hist(sample.means, cex.axis = 1.5, cex.lab = 1.5, main = "", freq = FALSE)
  c(mean(sample.means),
    var(sample.means), 
    sd(sample.means), 
    pi, 
    pi*(1-pi)/n, 
    sqrt(pi*(1-pi)/n))
}
dist.mean(8, 0.5)
[1] 0.50172500 0.03079385 0.17548178 0.50000000 0.03125000 0.17677670
dist.mean(32, 0.5)
[1] 0.50125625 0.00771267 0.08782181 0.50000000 0.00781250 0.08838835
dist.mean(128, 0.5)
[1] 0.499873438 0.001908915 0.043691133 0.500000000 0.001953125 0.044194174
dist.mean(512, 0.5)
[1] 0.5002300781 0.0004890104 0.0221135795 0.5000000000 0.0004882812
[6] 0.0220970869
dist.mean(2048, 0.5)
[1] 0.4999568848 0.0001207809 0.0109900359 0.5000000000 0.0001220703
[6] 0.0110485435

Beyond the binomial case?

  1. Think about other types of situations which do not fit the binomial case. They should be relatively easy to find.

    • Recall that a binomial random variable counts the number of “successes”.
    • Each trial results in “success” or “failure”.
    • The trials are independent.
    • The probability of success, \(\pi\), is constant across the trials.
    • There is a fixed number of trials, \(n\).
  1. Explore first what happens when \(n\to\infty\).

    • We already have seen what happens to a binomial random variable when \(n\to\infty\), holding everything else constant?
    • But we are going to dig further.
    • What happens to the shape, center, spread of the binomial distribution?
  1. There is a need to standardize, i.e. consider \[Z=\frac{X-\mathbb{E}\left(X\right)}{\mathsf{SD}\left(X\right)}=\frac{X-n\pi}{\sqrt{n\pi\left(1-\pi\right)}}\]

  2. This is the same as writing: \[Z=\frac{\widehat{\pi}-\pi}{\sqrt{\dfrac{\pi\left(1-\pi\right)}{n}}}\]

  1. What happens to the probability distribution of \(Z\) as \(n\to\infty\)?

  2. You will now see a Monte Carlo simulation of the behavior of \(Z\) as \(n\to\infty\).

    • Must generate simulated data via rbinom().
    • Must think about the underlying parameters to generate simulated data.
  3. Settings:

    • \(n\) could be 8, 32, 128, …
    • \(\pi=0.5\)
plot.dist.mean <- function(n, pi)
{
  a <- rbinom(10^4, n, pi)
  st.sample.means <- (a/n-pi)/(sqrt(pi*(1-pi)/n))
  par(mfrow = c(1, 2))
  hist(st.sample.means, cex.axis = 1.5, cex.lab = 1.5, main = "", xlim = c(-4,4), freq = FALSE)
  plot(ecdf(st.sample.means))
  curve(pnorm, add=TRUE, col="red")
}
plot.dist.mean(8, 0.5)
plot.dist.mean(32, 0.5)
plot.dist.mean(128, 0.5)
plot.dist.mean(512, 0.5)
plot.dist.mean(2048, 0.5)

What is the red curve?

  1. The red curve is actually the cumulative distribution function (cdf) of a random variable \(Z\) that has a standard normal distribution.

  2. The standard normal cdf is usually written at \(\Phi(z)\) has the following form: \[\Phi(z)=\mathbb{P}\left(Z\leq z\right)=\int_{-\infty}^z \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}x^2\right)\, dx.\]

  1. You will not be able to calculate the integral by hand. That is why standard normal tables and pnorm() are available.

  2. This is probably one of the more famous continuous distributions.

  3. What matters is that the cdf of a standard normal random variable is a good approximation of the cdf of the sampling distribution of the standardized sample mean. This is the central limit theorem in action.

Digression: Continuous random variables

In the slides, pay attention to whether the word “discrete” was specified or not. Differences compared to discrete random variables:

  1. Probability density functions (pdf) versus probability mass functions (pmf)

    1. Heights of those functions differ in meaning: for the pdf it would be density, for the pmf it would be probability
    2. Heights of a pdf can exceed 1.
  1. In terms of calculations, we have for a continuous random variable \(X\):

    • \(\mathbb{P}\left(X=x\right)=0\)
    • \(\displaystyle\mathbb{P}\left(a\leq X\leq b\right)=\mathbb{P}\left(a< X< b\right)=\int_a^b f_X(x)\,dx\)
    • \(\displaystyle\mathbb{E}\left(X\right)=\int_{-\infty}^\infty xf_X(x)\,dx\)

Digression: \(N(0,1)\)

The standard normal distribution has some nice properties:

  • Its pdf is symmetric about zero. The curve has a peak at 0.
  • One unit to the left and the right of zero: inflection points, PDF curve changes concavity
  • The area between -1 to 1 is roughly 68%. The area between -2 to 2 is roughly 95%. The area between -3 to 3 is roughly 99.7%.
  • Typically, probabilities from the standard normal are obtained from a table or from software.
par(mfrow=c(1,2)) # only for design purposes
plot(dnorm, from = -4, to = 4)
plot(pnorm, from = -4, to = 4)

What does the red curve buy us?

  1. We can now calculate probabilities involving the sample mean by directly using a large-sample approximation.

  2. To be more specific, suppose we want to know \(\mathbb{P}\left(X\geq 9\right)\) for \(X\sim \mathsf{Bin}\left(10,0.5\right)\) for Harley’s case.

  3. But you do not know simulation and you do not know anything about how to use the binomial distribution directly.

  1. Try applying the central limit theorem. How is \(X\) related to a standardized sample mean?

  2. Let \(X\sim \mathsf{Bin}\left(10,0.5\right)\).

    • Since the target is \(\mathbb{P}\left(X\geq 9\right)\), we can express this as \(1-\mathbb{P}\left(X \leq 8\right)\).
    • Whenever \(X\leq 8\), then \(Z\leq (0.8-0.5)/\sqrt{(0.5)(0.5)/10}\).
  1. Putting everything together, we have \[\begin{eqnarray*}\mathbb{P}\left(X\geq 9\right) &=& 1-\mathbb{P}\left(X \leq 8\right) \\ &=& 1-\mathbb{P}\left(Z \leq 1.90\right) \\ &\approx & 1- 0.97 = 0.03\end{eqnarray*}\]
1-pnorm((0.8-0.5)/sqrt(0.5*0.5/10))
[1] 0.02888979
  1. This is different compared to the what we got from simulation and from the exact calculations from the binomial.
  1. Does this mean the central limit theorem is useless?

    • Not really, you need \(n\to\infty\).
    • Here \(n=10\) only.
  2. One trick you can use for the binomial case is the continuity correction:

\[\begin{eqnarray*}\mathbb{P}\left(X\geq 9\right) &=& 1-\mathbb{P}\left(X \leq 8.5\right) \\ &=& 1-\mathbb{P}\left(Z \leq 2.21\right) \\ &\approx & 1- 0.99 = 0.01\end{eqnarray*}\]

1-pnorm((0.85-0.5)/sqrt(0.5*0.5/10))
[1] 0.01342835

The impact of a plug-in for \(\pi\)

  1. Recall the importance of the following random variable: \[Z=\frac{\widehat{\pi}-\pi}{\sqrt{\dfrac{\pi\left(1-\pi\right)}{n}}}\]

  2. Consider the alternative where \[Z=\frac{\widehat{\pi}-\pi}{\sqrt{\dfrac{\widehat{\pi}\left(1-\widehat{\pi}\right)}{n}}}=\frac{\widehat{\pi}-\pi}{\widehat{\mathsf{SE}}}\]

  1. How will the central limit theorem be affected?

  2. Recall the code to demonstrate the central limit theorem:

plot.dist.mean <- function(n, pi)
{
  a <- rbinom(10^4, n, pi)
  st.sample.means <- (a/n-pi)/(sqrt((a/n)*(1-a/n)/n))
  par(mfrow = c(1, 2))
  hist(st.sample.means, cex.axis = 1.5, cex.lab = 1.5, main = "", xlim = c(-4,4), freq = FALSE)
  plot(ecdf(st.sample.means))
  curve(pnorm, add=TRUE, col="red")
}
  1. Observe what changed in the code.
plot.dist.mean(32, 0.5)
plot.dist.mean(128, 0.5)
plot.dist.mean(512, 0.5)
plot.dist.mean(2048, 0.5)
  1. Intuitively, \(\widehat{\pi}\) becomes practically \(\pi\) when \(n\to\infty\).
  2. As a result, the random variable \[Z=\frac{\widehat{\pi}-\pi}{\sqrt{\dfrac{\widehat{\pi}\left(1-\widehat{\pi}\right)}{n}}}=\frac{\widehat{\pi}-\pi}{\widehat{\mathsf{SE}}}\] can be used for hypothesis testing as well.

Can we obtain confidence intervals for \(\pi\) via the central limit theorem?

  1. Yes and these intervals become easier to compute.

  2. Goal is the same as before. We want another way of obtaining a \(100(1-\alpha)%\) confidence interval for \(\pi\).

  3. The approach is the same as before, but we let \(x_{\mathsf{obs}}\) be the observed value of \(X\sim \mathsf{Bin}\left(n,\pi\right)\).

    • Adjust the value of \(\pi\) so that \(\mathbb{P}\left(X\geq x_{\mathsf{obs}}\right) \approx \alpha/2\).
    • Adjust the value of \(\pi\) so that \(\mathbb{P}\left(X\leq x_{\mathsf{obs}}\right) \approx \alpha/2\).
  1. Start from \(\mathbb{P}\left(X\leq x_{\mathsf{obs}}\right)\). Express this in terms of \(Z\), i.e. \[\begin{eqnarray*} \mathbb{P}\left(X\leq x_{\mathsf{obs}}\right) = \mathbb{P}\left(Z\leq \frac{x_{\mathsf{obs}}/n-\pi}{\widehat{\mathsf{SE}}}\right) \end{eqnarray*}\]

  2. The target is to have \[\begin{eqnarray*} \mathbb{P}\left(X\leq x_{\mathsf{obs}}\right) = \mathbb{P}\left(Z\leq \frac{x_{\mathsf{obs}}/n-\pi}{\widehat{\mathsf{SE}}}\right) \approx \alpha/2 \end{eqnarray*}\] and solve for \(\pi\).

  1. Solving for \(\pi\) requires solving a simpler problem: Find a number \(z_{\alpha/2}\) such that \[\mathbb{P}\left(Z\leq z_{\alpha/2}\right) \approx \alpha/2\]

  2. This can be obtained easily from a standard normal table or from qnorm().

  3. Once you found \(z_{\alpha/2}\), solve for the upper limit of \(\pi\) in \[\frac{x_{\mathsf{obs}}/n-\pi}{\widehat{\mathsf{SE}}}=z_{\alpha/2}\]

  1. Similarly, \[\begin{eqnarray*} \mathbb{P}\left(X\geq x_{\mathsf{obs}}\right) &=& \mathbb{P}\left(Z\geq \frac{x_{\mathsf{obs}}/n-\pi}{\widehat{\mathsf{SE}}}\right) \\ &=& 1-\mathbb{P}\left(Z\leq \frac{x_{\mathsf{obs}}/n-\pi}{\widehat{\mathsf{SE}}}\right) \approx \alpha/2 \end{eqnarray*}\] and solve for \(\pi\).

  2. Solving for \(\pi\) here requires solving a simpler problem: Find a number \(z_{1-\alpha/2}\) such that \[\mathbb{P}\left(Z\leq z_{1-\alpha/2}\right) \approx 1-\alpha/2\]

  1. Once you found \(z_{1-\alpha/2}\), solve for the lower limit of \(\pi\) in \[\frac{x_{\mathsf{obs}}/n-\pi}{\widehat{\mathsf{SE}}}=z_{1-\alpha/2}\]

  2. Because of the symmetry of the standard normal distribution, \(z_{1-\alpha/2}=-z_{\alpha/2}\).

  3. Therefore, \[x_{\mathsf{obs}}/n +z_{\alpha/2}\widehat{\mathsf{SE}}\leq \pi \leq x_{\mathsf{obs}}/n -z_{\alpha/2}\widehat{\mathsf{SE}}\]

  1. We can definitely treat the lower and upper limits of a confidence interval as random variables so that we have a procedure for constructing confidence intervals with stated performance guaranetees:

\[\widehat{\pi} +z_{\alpha/2}\widehat{\mathsf{SE}}\leq \pi \leq \widehat{\pi} -z_{\alpha/2}\widehat{\mathsf{SE}}\]

  1. This interval is a normal-based large-sample (or asymptotic) 95% confidence interval for \(\pi\).

  2. What we did in earlier slides does not require \(n\to\infty\).

Moving forward

  1. The central limit theorem applies to standardized sample means, even if the data do not necessarily consist of 0’s and 1’s.

  2. Let \(X_1,\ldots, X_n\) be random variables and that \(\overline{X}\) is the sample mean.

  3. \(i\) is the index to represent which draw from a population or process.

  1. The main restrictions are (mirroring the binomial case):

    • The draws are independent.
    • The draws come from the same distribution, even if we do not know the actual distribution.
    • The target parameter of interest is the common expected value \(\mu=\mathbb{E}\left(X_i\right)\).
    • The variance of each random variable \(\sigma^2=\mathsf{Var}\left(X_i\right)\) is finite.
  1. In the binomial case:

    • \(X_1, X_2, \ldots, X_n\) independent
    • \(X_1, X_2, \ldots, X_n\) have the same distribution
    • \(\mathbb{E}\left(X_i\right)=\pi\)
    • \(\mathsf{Var}\left(X_i\right)=\pi\left(1-\pi\right)\)
  2. What happens if we want to use simulation-based approaches, outside of the binomial case?

  1. Every time you apply statistics to any research question, you always think about the following:

    • What is the parameter?
    • How do I estimate the parameter from the data?
    • Can I express uncertainty resulting from sampling variation? Can I compute this from the data?
    • Can I produce a confidence interval for the parameter?
    • Can I test hypothesis about the parameter?
  1. Think of every thing you do with the data as a procedure.

    • Does the procedure have desirable properties or performance guarantees?
    • What conditions are needed for the pocedure to attain these desirable properties or performance guarantees?
  2. Once you think this way, you can apply statistics to almost any situation where there is sampling variation to grapple with.