Distribution of random sample of confidence intervals containing the true mean Given 50 random samples, each of size 25, from a normal distribution with mean 20 and standard deviation 5. From each of the 50 samples, you can find a 90% confidence interval for the mean. Let Y be a random variable indicating how many of the confidence intervals containing the mean. What is the distribution of Y?
My friend and I are having a heated discussion about this. I think it is binomially distributed with 50 observations and probability of success 90%. My friend claims it is normally distributed with mean 45 and standard deviation 1.
I think this because Y has a fixed number of observations (50), each observation has two outcomes, and their probabilities are fixed (probability of containing the mean is 90%).
 A: Your reasoning is good.  The purpose of this answer is to show how you can check it out.  It's a very useful skill, because we all make mistakes--and this allows us to find them in private before they become embarrassing :-)
The process described below is the important thing to follow, because it can be used over and over again to address all kinds of statistical questions, both simple and sophisticated.  I am using R because it is free and popular, but the same approach will work in other statistical computing platforms, such as Python, Mathematica, Matlab, or whatever.
The heart of the question concerns a confidence interval for a sample from a Normal distribution.  So, we begin with code to compute a CI of a set of numbers z and check whether it covers the mean mu with confidence level 1-alpha:
covers <- function(z, mu=0, alpha=0.10) {
  n <- length(z)
  abs(mean(z) - mu) <= qt(1 - alpha/2, n-1) * sqrt(var(z) / n)
}

This computation should look familiar, involving the mean, variance, and length of z along with a comparison to a critical Student t-value computed with qt.  The comparison results in a value of 1 when true and 0 otherwise: this is the output of covers.
To simulate one experiment we need to (1) generate random samples z; (2) do this 50 times over; and (3) count how many of those 50 samples have CIs that cover the mean.  As in much of mathematics and computer science, the correct expression is developed from the inside out, so the following line executes these operations in the order (3)-(2)-(1) using sum (which, by summing the zeros and ones returned by covers, is precisely the count associated with the random variable $Y$), apply, and rnorm, respectively:
n <- 50
k <- 25
mu <- 20
sum(apply(matrix(rnorm(k*n, mu, 5), nrow=k), 2, covers, mu=mu))

To make it a little more readable and a lot more flexible I have placed the numbers 50, 25, and 20 in variables.  Changing the values of these variables lets you play with variations of the problem.
Running a simulation calls for repeating this operation over and over and collecting the results.  Do this by wrapping the preceding within a call to replicate:
sim <- replicate(1000, {
  sum(apply(matrix(rnorm(k*n, mu, 5), nrow=k), 2, covers, mu=mu))
})

The argument 1000 specifies how many independent runs to perform.  Because I have made no effort at computational efficiency, this will take a few seconds for a thousand iterations.  But a thousand will be more than enough.  All these results are collected into the array named sim for further analysis.
Finally, we should plot the results and perhaps perform a few tests, as if such formalities really mattered:
hist(sim, freq=FALSE, breaks=min(sim):(n+1)-1/2, ylim=c(0, dnorm(0.9*n, 0.9*n, 1)),xlab="Y")
curve(dnorm(x, 0.9*n, 1), col="Red", add=TRUE)
points(min(sim):n, dbinom(min(sim):n, n, 0.9), pch=19, col="Blue")

Over a histogram (via hist) I have superimposed the putative Normal density function (via curve) as well as the (correct) Binomial distribution (via points).  Clearly the Normal distribution isn't even a good approximation to the results: although its peak is in the right place, its standard deviation is far too small.

A $\chi^2$ test amply supports the visual impression that the binomial distribution (blue dots) is the true underlying distribution:
chisq.test(tabulate(sim+1, n+1), p=dbinom(0:n, n, 0.9), simulate.p.value=TRUE)

The output is

X-squared = 10.0457, df = NA, p-value = 0.7441

The large p-value says the simulation results cannot be distinguished from 1000 independent draws from a Binomial$(50, 0.9)$ distribution.  (Your output will differ slightly because the random values in your simulation will differ from mine.)
