\documentclass[12pt]{article} \usepackage{Math203} \usepackage{amsmath,amssymb} \usepackage{hyperref} \usepackage{color} \usepackage{xcolor} \usepackage{enumitem} \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{The Geometric and Poisson Distributions} \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 \textbf{Geometric Distribution:} The Geometric distribution is a probability distribution for a random variable derived from sequences of \textit{binary trials}. As for the binomial distribution, wach trial has `success probability' $p$, that is, the probability that the outcome is 1 is $p$, and the probability that the outcome is 0 is $q=1-p$, \bigskip \textbf{Experimental set up:} We consider the experiment involving a sequence of identical binary trials with success probability $p$ whose outcomes are independent. The experiment terminates when a success is observed. For example, we might have five experimental results of results \[ 00001\\ 001\\ 1\\ 0001\\ 00000001 \] corresponding to the experiment terminating after 5, 3, 1, 4 and 8 trials respectively. \bigskip \textbf{Random variable:} To define the random variable, $Y$, we record the \textbf{number of trials} until the first success is obtained. If $0 < p < 1$, the set of possible values $Y$ can take is \[ \calY = \{1,2,3,\ldots,\}. \] with no upper limit on the possible values. All of these values will carry some positive probability. \bigskip \textbf{Probability distribution:} This definition for $Y$ allows us to compute the probability that should be attached to each of the possible values that $Y$ can take. We consider events of the form $(Y=y)$ for numerical values $y=1,2,3,\ldots,$, and observe that \[ (Y = y) \quad \text{means} \quad ``\text{first successes on trial $y$}" \] In the examples above, the probability associated with each sequence is \[ \begin{array}{ll} 00001 & qqqqp\\ 001 & qqp\\ 1 & p\\ 0001 & qqqp\\ 00000001 & qqqqqqqp. \end{array} \] Therefore we have that the probability mass function (pmf) for the Geometric distribution is \[ p(y) = P(Y=y) = \left\{ \begin{array}{ll} q^{y-1} p & y \in \{1,2,3,\ldots,\}\\[6pt] 0 & \text{otherwise} \end{array} \right. \] We write \[ Y \sim Geometric(p). \] \pagebreak \textbf{Numerical calculations:} Numerical calculations for the Geometric distribution can be performed in \texttt{R}. For $p=0.4$, <>= p<-0.4 q<-1-p y<-3 prob<-q^{y-1}*p print(paste('P(Y=3) =',prob)) @ The calculation can also be performed using the specially written function \texttt{dgeom} <>= print(paste('P(Y=3) =',dgeom(y-1,p))) @ Note that you need to plug $y-1$ (not $y$) into the function, due to the way the \texttt{dgeom} function is written. We can compute the complete list of values for the pmf evaluated at $y=1,2,\ldots,10$. <>= y<-1:10; probs<-dgeom(y-1,p) print(cbind(y,probs)) @ We can plot the pmf as a function of $y$. <>= par(mar=c(4,4,2,2)) plot(y,probs,pch=19,ylab='p(y)',xlim=range(0,10),ylim=range(0,0.6),xaxt='n') axis(1,at=0:10) @ \medskip \textbf{Connection to samples:} We can simulate data mimicking the experiment. For example, we can generate a matrix of sequences of binary trials, and then compute $Y$ for each sequence; here for $N=5$ sequences <>= set.seed(3) N<-5 Yseq<-list() for(i in 1:5){ success<-FALSE bseq<-as.character() while(!success){ X<-rbinom(1,1,p) bseq<-paste(bseq,X,sep='') if(X == 1) success<-TRUE } Yseq[[i]]<-bseq print(bseq) } Y<-as.numeric(lapply(Yseq,nchar)) Y @ Hence the simulated $Y$ values are \Sexpr{Y}. We can also simulate $Y$ directly using \texttt{rgeom}. Here are $N=10$ values <>= N<-10;rgeom(N,p)+1 @ Suppose we simulate a very large sample $N=1000000$, and draw a density histogram <>= N<-1000000 Y<-rgeom(N,p)+1 par(mar=c(4,4,0,0)) hist(Y,breaks=0:max(Y),right=FALSE,include.lowest=TRUE,main='', xaxt='n',freq=FALSE,ylim=range(0,0.6)) axis(1,at=c(0:max(Y))) box() y<-1:max(Y) probs<-dgeom(y-1,p) points(y,probs,pch=19,col='red') legend(6,0.3,c('p(y)'),pch=19,col='red') @ The pmf values are superimposed in red. The pmf can therefore be interpreted as the limiting values of the density histogram heights in a very large sample. \bigskip In the sample of \Sexpr{N} from the $Geometric(\Sexpr{p})$ distribution, we have the following counts of the different possible values. <>= table(Y) @ The relative frequencies are <>= table(Y)/N @ which closely match the pmf values <>= round(cbind(table(Y)/N,probs),6) @ \textbf{Expectation and Variance:} If $Y\sim Geometric(p)$, then the (population) mean (or expectation) and variance are \[ E(Y) = \frac{1}{p} \qquad Var(Y) = \frac{(1-p)}{p^2} = \frac{q}{p^2} \] <>= c(1/p,q/p^2) #Expectation, Variance @ That is, the expectation is $\Sexpr{1/p}$ and the variance is $\Sexpr{q/p^2}$. These values describe aspects of the distribution, and the function $p(y)$. In the very large sample of $N=100000$ values, we have sample mean and sample variance which are very close to the population expectation and variance <>= mean(Y) var(Y) @ \pagebreak \textbf{The Poisson distribution:} The Poisson distribution arises as a limiting case of the Binomial distribution when $n$ is large and $p$ is small. \bigskip \textbf{Experimental set up:} We consider the experiment involving a sequence of $n$ identical binary trials with success probability $p$ whose outcomes are independent. Suppose that $n$ is large and $p$ is small, and let $\lambda = np$. We consider the limiting case where $n \longrightarrow \infty$ and $p \longrightarrow 0$ but where $\lambda$ is held fixed: that is, we fix $\lambda$ to be some positive value, and let \[ p = \frac{\lambda}{n}. \] \bigskip \textbf{Random variable:} We record the \textbf{total number} of successes in the sequence. The set of possible values $Y$ can take is \[ \calY = \{0,1,2,3,\ldots\}. \] where $n$ is so large that we might as well consider it infinite. \bigskip \textbf{Probability distribution:} We consider events of the form $(Y=y)$ for numerical values $y=0,1,2,\ldots,$, and observe that \[ (Y = y) \quad \text{means} \quad ``\text{$y$ successes}" \] It can be shown that \[ p(y) = P(Y=y) = e^{-\lambda} \frac{\lambda^y}{y!} \qquad y \in \{0,1,2,\ldots,\} \] with $p(y)=0$ otherwise. We write \[ Y \sim Poisson(\lambda). \] \bigskip \textbf{Numerical calculations:} Numerical calculations for the Poisson distribution can be performed in \texttt{R}. We use the settings $\lambda=2$ and $n=100000$ <>= lambda<-2 n<-100000 p<-lambda/n y<-3 prob<-choose(n,y)*(p^y)*((1-p)^(n-y)) #choose(.,.) computes C_n^r. print(paste('P(Y=3) =',round(prob,6))) @ The calculation can also be performed using \texttt{dbinom} <>= print(paste('P(Y=3) =',round(dbinom(y,n,p),6))) @ In fact, there is a specially written function \texttt{dpois} to deal with this case <>= print(paste('P(Y=3) =',round(dpois(y,lambda),6))) @ We can compute the complete list of values for the pmf evaluated at $y=0,1,2,\ldots,10$. <>= y<-0:10; probs<-dbinom(y,n,p); probsPois<-dpois(y,lambda) print(cbind(y,round(probs,6), round(probsPois,6))) @ We can plot the pmf as a function of $y$. <>= par(mar=c(4,4,2,2)) plot(y,probs,pch=19,ylab='p(y)',ylim=range(0,0.6),xaxt='n');axis(1,at=0:10) points(y,probsPois,pch=3,col='red') legend(8,0.5,c('Binomial(n,p)',expression(Poisson(lambda))), pch=c(19,3),col=c('black','red')) @ \medskip \textbf{Connection to samples:} We can also $Y$ directly using \texttt{rpois}. Here are $N=10$ values <>= N<-10;rpois(N,lambda) @ Suppose we simulate a very large sample $N=10000$, and draw a density histogram <>= N<-10000 Y<-rpois(N,lambda) par(mar=c(4,4,0,0)) hist(Y,breaks=0:10,right=FALSE,include.lowest=TRUE,main='', xaxt='n',freq=FALSE,ylim=range(0,0.6)) axis(1,at=c(0:10)) box() points(y,probs,pch=19,col='red') legend(6,0.3,c('p(y)'),pch=19,col='red') @ The pmf values are superimposed in red. The pmf can again be interpreted as the limiting values of the density histogram heights in a very large sample. \bigskip In the sample of \Sexpr{N} from the $Poisson(\Sexpr{lambda})$ distribution, we have the following counts of the different possible values. <>= table(Y) @ The relative frequencies are <>= table(Y)/N @ which closely match the pmf values <>= round(cbind(table(Y)/N,probsPois),6) @ \textbf{Expectation and Variance:} If $Y\sim Poisson(\lambda)$, then the (population) mean (or expectation) and variance are \[ E(Y) = \lambda \qquad Var(Y) = \lambda \] That is, the expectation and variance are both $\Sexpr{lambda}$. In the very large sample of $N=100000$ values, we have sample mean and sample variance which are very close to the population expectation and variance <>= mean(Y) var(Y) @ \end{document}