\documentclass[12pt]{article} \usepackage{Math203} \usepackage{amsmath,amssymb} \usepackage{hyperref} \usepackage{color} \usepackage{xcolor} \usepackage{enumitem} \usepackage{longtable} \usepackage{booktabs} \pagestyle{plain} \renewcommand{\labelenumi}{\arabic{enumi}.} \renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}} \renewcommand{\labelenumiii}{\arabic{enumi}.\arabic{enumii}.\arabic{enumiii}} \setlength{\topmargin}{-0.75 in} \setlength{\textheight}{9.5 in} \parindent0in \newenvironment{packed_enum}{ \begin{enumerate} \setlength{\itemsep}{12pt} \setlength{\parskip}{0pt} \setlength{\parsep}{0pt} \setlength{\topsep}{0pt} }{\end{enumerate}} \renewcommand{\labelitemii}{$\blacktriangleright$} \def\ctitle{Math 203: Principles of Statistics 1} \def\ttitle{Fall 2024} \def\stitle{Using the Central Limit Theorem} \makeatletter \def\@maketitle{% \newpage \null \begin{center}% {\textsclarge{\ctitle} \par}% \medskip {\textsclarge{\ttitle} \par}% \medskip {\textsc{\stitle} \par} \end{center}% \par \vskip 1.5em} \makeatother \begin{document} <>= myrnd <- function(x,d) {formatC(x, format="f", digits=d)} library(shape) library(diagram) @ \maketitle The Central Limit Theorem can be used to approximate the distributions of sums and means of independent and identically distributed (i.i.d.) random variables. If $X_1,\ldots,X_n$ are i.i.d. (that is, they are a random sample) from a distribution with expectation $\mu$ and variance $\sigma^2$, then \[ S_n = \sum_{i=1}^n X_i \qquad \text{and} \qquad \Xbar = \frac{1}{n}\sum_{i=1}^n X_i \] are random variables with that are approximately Normally distributed. That is, \[ S_n \approxsim Normal(n\mu, n \sigma^2) \qquad \text{and} \qquad \Xbar \approxsim Normal(\mu, \sigma^2/n) \] Therefore, \[ P(S_n \leq s) \approx \Phi\left(\frac{s - n \mu}{\sqrt{n \sigma^2}} \right) \qquad P(\Xbar \leq x) \approx \Phi\left(\frac{x - \mu}{\sigma/\sqrt{n}} \right) \] where $\Phi(\cdot)$ is the standard Normal cdf. \bigskip \textbf{Example:} Suppose $X$ records the score when a die is rolled once. Then $X$ has a discrete Uniform distribution on the values $\{1,2,3,4,5,6\}$ and \[ P(X = x) = \frac{1}{6} \qquad x=1,2,3,4,5,6 \] with zero probability for any other value of $x$. We have that \[ E(X) = \sum_{x=1}^6 x p(x) = \sum_{x=1}^6 x \times \frac{1}{6} = \frac{21}{6} = \Sexpr{x<-1:6;mu<-sum(x)/6;sigsq<-sum(x^2/6)-(sum(x)/6)^2;myrnd(mu,1)} \] and \[ E(X^2) = \sum_{x=1}^6 x^2 p(x) = \sum_{x=1}^6 x^2 \times \frac{1}{6} = \frac{91}{6} = \Sexpr{myrnd(sum(x^2/6),4)} \] so that \[ Var(X) = E(X^2) - \{E(X)\}^2 = \frac{91}{6} - \left(\frac{21}{6}\right)^2 = \Sexpr{x<-1:6;mu<-sum(x)/6;sigsq<-sum(x^2/6)-(sum(x)/6)^2;myrnd(sigsq,4)}. \] Suppose now we roll the die $n=50$ times independently, with $X_1,\ldots,X_n$ recording the individual scores. Then the total score across all rolls is $S_n$, and the Central Limit Theorem says that \[ S_n \approxsim Normal(n\mu, n \sigma^2) = Normal(50 \times \Sexpr{myrnd(mu,1)}, 50 \times \Sexpr{myrnd(sigsq,4)}) = Normal(\Sexpr{myrnd(50*mu,1)}, \Sexpr{myrnd(50*sigsq,4)}) \] and for the sample mean score $\Xbar$ \[ \Xbar \approxsim Normal(\mu, \sigma^2/n) = Normal(\Sexpr{myrnd(mu,1)}, \Sexpr{myrnd(sigsq,4)}/50) = Normal(\Sexpr{myrnd(mu,1)}, \Sexpr{myrnd(sigsq/50,4)}) \] Then for example, using the Z-transform approach, we have \[ P(S_n > 190) \approx P\left(\frac{S_n - \Sexpr{50*mu}}{\sqrt{\Sexpr{myrnd(50*sigsq,4)}}} > \frac{190 - \Sexpr{50*mu}}{\sqrt{\Sexpr{myrnd(50*sigsq,4)}}} \right) = P(Z_n > \Sexpr{myrnd((190-50*mu)/sqrt(50*sigsq),4)}) \] and so \[ P(S_n > 190) = 1 - P(Z_n \leq \Sexpr{myrnd((190-50*mu)/sqrt(50*sigsq),4)}) = 1 - \Sexpr{myrnd(pnorm((190-50*mu)/sqrt(50*sigsq)),4)} = \Sexpr{myrnd(1-pnorm((190-50*mu)/sqrt(50*sigsq)),4)} \] We can test this in simulation, using $N=100000$ repeat runs of the experiment. The following histogram summarizes the observed values of $S_n$ over the repeat runs, with the computed Normal approximation. <>= set.seed(2) N<-100000;n<-50 par(mar=c(4,4,2,2)) Sn<-replicate(N,sum(sample(1:6,size=n,rep=T))) tty<-expression(paste('Density histogram for ', S[n], ' over 100000 repeats')) hist(Sn,br=seq(110,230,by=1),main=tty, freq=FALSE,xlab=expression(S[n]));box() xv<-seq(110,230,by=0.01);yv<-dnorm(xv,n*mu,sqrt(n*sigsq));lines(xv,yv,col='red') legend(110,0.30,c('Normal approximation'),col='red') @ In the repeat runs, the number of times $S_n$ exceeds 190 is \Sexpr{sum(Sn>190)}, so the estimated probability of exceeding 190 is \Sexpr{sum(Sn>190)/N}, which is fairly close to the claimed value \Sexpr{myrnd(1-pnorm((190-50*mu)/sqrt(50*sigsq)),4)}. \pagebreak In this problem, the sum random variable $S_n$ is \textbf{discrete}; there are $6^{50}$ possible sequences of rolls, and the resulting sum can take any whole number value between 60 and 300. If we consider the required probability, then it is evident that in the exact version \[ P(S_n > 190) = P(S_n \geq 191) = 1 - P(S_n \leq 190) \] If we use the Central Limit Theorem to compute the two equivalent versions, we get two slightly different answers: <>= n<-50 mu<-sum(x)/6 sigsq<-sum(x^2/6)-(sum(x)/6)^2 1-pnorm(190,n*mu,sqrt(n*sigsq)) 1-pnorm(191,n*mu,sqrt(n*sigsq)) @ In such cases, a \textit{continuity correction} is sometimes applied, and we compute instead \[ P(S_n > 190.5) \] adding on 0.5 to the $y$ value <>= 1-pnorm(190.5,n*mu,sqrt(n*sigsq)) @ \textbf{Example:} Suppose that $Y \sim Binomial(n,p)$; recall that $Y$ can be viewed as the total number of successes on $n$ independent binary trials with success probability $p$. We may therefore write \[ Y = \sum_{i=1}^n X_i \] where $X_i$ is the result of the $i$th trial. We may then apply the Central Limit Theorem to $Y$ (as it is the sum of independent and identically distributed random variables) and say that \[ Y \approxsim Normal(np, npq) \] as $E(X_i) = p$ and $Var(X_i) = p (1-p) = pq$. For a numerical illustration, suppose $n=500$ and $p=0.4$. We can compute exactly that $P(Y > 210) = \Sexpr{1-pbinom(210,500,0.4)}$. <>= n<-500 p<-0.4 q<-1-p y<-210 1-pbinom(y,n,p) @ For a Central Limit Theorem approximation \[ P(Y > 210) = P\left(\frac{Y - np}{\sqrt{npq}} > \frac{210 - np}{\sqrt{npq}} \right) = P\left(\frac{Y - \Sexpr{n}\times\Sexpr{p}}{\sqrt{\Sexpr{n}\times\Sexpr{p}\times\Sexpr{q}}} > \frac{210 - \Sexpr{n}\times\Sexpr{p}}{\sqrt{\Sexpr{n}\times\Sexpr{p}\times\Sexpr{q}}} \right) \] or \[ P(Y > 210) = P(Z > \Sexpr{(y-n*p)/sqrt(n*p*q)}) = 1-P(Z \leq \Sexpr{(y-n*p)/sqrt(n*p*q)}) \approx 1-\Sexpr{myrnd(pnorm((y-n*p)/sqrt(n*p*q)),4)} = \Sexpr{myrnd(1-pnorm((y-n*p)/sqrt(n*p*q)),4)}. \] This value is not a very close approximation. Using the continuity correction \[ P(Y > 210) \approx P\left(\frac{Y - np}{\sqrt{npq}} > \frac{210+0.5 - np}{\sqrt{npq}} \right) = P\left(\frac{Y - \Sexpr{n}\times\Sexpr{p}}{\sqrt{\Sexpr{n}\times\Sexpr{p}\times\Sexpr{q}}} > \frac{210.5 - \Sexpr{n}\times\Sexpr{p}}{\sqrt{\Sexpr{n}\times\Sexpr{p}\times\Sexpr{q}}} \right) \] or \[ P(Y > 210) \approx P(Z > \Sexpr{(y+0.5-n*p)/sqrt(n*p*q)}) = 1-P(Z \leq \Sexpr{(y+0.5-n*p)/sqrt(n*p*q)}) \approx 1-\Sexpr{myrnd(pnorm((y+0.5-n*p)/sqrt(n*p*q)),4)} = \Sexpr{myrnd(1-pnorm((y+0.5-n*p)/sqrt(n*p*q)),4)}. \] This value is much closer to the exact value. Similarly, suppose we wish to compute \[ P(Y \leq 195) \] which in the exact version equals $\Sexpr{pbinom(195,n,p)}$. <>= pnorm(195,n*p,sqrt(n*p*q)) #Without continuity correction pnorm(195.5,n*p,sqrt(n*p*q)) #With continuity correction @ The logic behind the continuity correction stems from the following diagram, which plots the binomial probabilities and the Normal approximation: <>= yv<-seq(0,n,by=1) pyv<-dbinom(yv,n,p) par(mar=c(4,4,2,2)) plot(yv,pyv,pch=19,cex=0.5, xlab='y', ylab='p(y)') title('Binomial(500,0.4) pmf') xv<-seq(0,n,by=0.001) fxv<-dnorm(xv,n*p,sqrt(n*p*q)) lines(xv,fxv,col='red') legend(0,0.035,c('Normal approximation'),col='red',lty=1) @ Overall, the Normal approximation to the Binomial does a good job. However, if we zoom into the region near the target value $y=190$ we see that the Normal approximation slightly underestimates the discrete probabilities <>= yv<-seq(190,200,by=1) pyv<-dbinom(yv,n,p) par(mar=c(4,4,3,2)) plot(yv,pyv,pch=19,cex=0.75, xlab='y', ylab='p(y)',xaxt='n') axis(1,at=yv) title('Binomial(500,0.4) pmf') xv<-seq(190,200,by=0.001) fxv<-dnorm(xv,n*p,sqrt(n*p*q)) lines(xv,fxv,col='red') legend(190,0.035,c('Normal approximation'),col='red',lty=1) @ \pagebreak If we look at the cumulative probabilities, the picture is clearer: <>= yv<-seq(190,200,by=1) pyv<-pbinom(yv,n,p) par(mar=c(4,4,3,2)) plot(yv,pyv,pch=19,cex=0.75,xlab='y', ylab='p(y)',xaxt='n') axis(1,at=yv) title('Binomial(500,0.4) cdf') xv<-seq(190,200,by=0.001) fxv<-pnorm(xv,n*p,sqrt(n*p*q)) lines(xv,fxv,col='red') legend(190,0.45,c('Normal cdf approximation'),col='red',lty=1) @ Therefore, to approximate the required cumulative probability, we need to increase the $y$ value by a small amount when using the Normal approximation. \bigskip As a final comparison, if we compute the difference \[ F(y)-\Phi \left( \frac{y-np}{\sqrt{npq}} \right) \] between the $Binomial(n,p)$ cdf and the $Normal(np,npq)$ cdf, we see that the difference is always positive, and peaks near $\mu = np$, whereas if we use the continuity correction \[ F(y)-\Phi \left( \frac{y+0.5-np}{\sqrt{npq}} \right) \] then the difference is smaller in magnitude, and hence the approximation is better. We plot this below for $y$ in the range 150 to 250. <>= yv<-seq(150,250,by=1) pyv<-pbinom(yv,n,p) fxv<-pnorm(yv,n*p,sqrt(n*p*q)) fxv.cc<-pnorm(yv+0.5,n*p,sqrt(n*p*q)) par(mar=c(4,4,3,2)) plot(yv,pyv-fxv,pch=19,cex=0.5,xlab='y', ylab='Difference in cdfs' ,ylim=range(-0.005,0.02)) points(yv,pyv-fxv.cc,col='green',pch=19,cex=0.5) title('Difference in cdfs (p=0.4)') legend(150,0.02,c('Without continuity correction','With continuity correction'),pch=19,col=c('black','green')) @ The same correction works well for other choices of $p$, provided the value is not very close to zero or one. <>= p<-0.75 q<-1-p yv<-seq(300,450,by=1) pyv<-pbinom(yv,n,p) fxv<-pnorm(yv,n*p,sqrt(n*p*q)) fxv.cc<-pnorm(yv+0.5,n*p,sqrt(n*p*q)) par(mar=c(4,4,3,2)) plot(yv,pyv-fxv,pch=19,cex=0.5,xlab='y', ylab='Difference in cdfs',ylim=range(-0.005,0.02)) points(yv,pyv-fxv.cc,col='green',pch=19,cex=0.5) title('Difference in cdfs (p=0.75)') legend(300,0.02,c('Without continuity correction','With continuity correction'),pch=19,col=c('black','green')) abline(v=n*p,h=0,lty=2) @ \pagebreak \textbf{Example:} Suppose $Y \sim Binomial(n,p)$ with $n=20$ and $p=0.6$. Here $n$ is much smaller, and we can investigate the continuity correction in more detail. If we draw the pmf of $Y$ as a bar plot, we obtain the following: <>= n<-20 p<-0.6; q<-1-p y<-0:n; py<-dbinom(y,n,p) mu<-n*p; sigsq<-n*p*q par(mar=c(4,4,3,2)) barplot(py ~ y,col='gray', xlab="y", ylab="p(y)", space=0, main="pmf for Binomial(20,0.6)",ylim=range(0,0.25)); box() xv<-seq(0,20,by=0.01); yv<-dnorm(xv,mu,sqrt(sigsq)) lines(xv+0.5,yv, col='red') abline(v=10.5,lty=2);legend(0,0.20,c('Normal approximation'),col='red',lty=1) @ The bar is plotted with bars centered at the labelled $y$ values. The Normal approximation here would be $Normal(np, npq)$, or $Normal(\Sexpr{n*p},\Sexpr{n*p*(1-p)})$, which is overlaid on the bar plot in red. We can see from the plot that there is always part of the each bar that is above the normal approximation curve; to compute $P(Y \leq 10)$, we use the Normal approximation up to $y$ value 10.5. <>= pbinom(10,n,p) #Exact binomial probability pnorm(10,n*p,sqrt(n*p*q)) #Normal approx. without continuity correction pnorm(10.5,n*p,sqrt(n*p*q)) #Normal approx. with continuity correction @ \end{document}