\documentclass[notitlepage]{article} \usepackage{../Math323} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{amsfonts,amsmath} \setlength{\parindent}{0pt} \def\E{\Expect} \lstloadlanguages{R} \definecolor{keywordcolor}{rgb}{0,0.6,0.6} \definecolor{delimcolor}{rgb}{0.461,0.039,0.102} \definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432} \lstdefinestyle{Rsettings}{ basicstyle=\ttfamily, breaklines=true, showstringspaces=false, keywords={if, else, function, theFunction, tmp}, % Write as many keywords otherkeywords={}, commentstyle=\itshape\color{Rcommentcolor}, keywordstyle=\color{keywordcolor}, moredelim=[s][\color{delimcolor}]{"}{"}, } \lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) @ \begin{center} {\textsclarge{Discrete Distributions: Examples in \texttt{R}}} \end{center} \begin{enumerate}[1.] \item \textbf{The Binomial distribution.} The binomial pmf $p(y)$ takes the form \[ p(y) = \dbinom{n}{y} p^y (1-p)^{n-y} \qquad y = 0,1,2,\ldots,n \] with $p(y) = 0$ for other values of $y$. This pmf can be computed using the \texttt{dbinom} function in \texttt{R}. <>= n<-20;p<-0.2;y<-0:n #Set the parameters p.y<-dbinom(y,n,p) #Compute the pmf par(mar=c(4,4,0,2)) #Set up the plot plot(y,p.y,pch=19,ylab='p(y)',cex=0.75) #Produce the plot print(cbind(y,p.y)) #Print the numerical values of y and p(y) @ We can use the function \texttt{pbinom} to compute the cdf $F(y)$ \[ F(y) = P(Y \leq y) = \sum_{j=0}^{\lfloor y \rfloor} p(j) = \sum_{j=0}^{\lfloor y \rfloor} \dbinom{n}{j} p^j (1-p)^{n-j} \qquad y > 0 \] where $\lfloor y \rfloor$ is the largest integer no greater than $y$. <>= F.y<-pbinom(y,n,p) #Compute the cdf par(mar=c(4,4,0,2)) plot(stepfun(y,c(0,F.y)),xlab='y',ylab='F(y)',pch=19,cex=0.5,vertical=FALSE,main='') @ We can compute the expectation and variance by using the \texttt{sum} function <>= Expectation<-sum(y*p.y) #Compute the expectation Variance<-sum((y-Expectation)^2*p.y) #Compute the variance c(Expectation,Variance) #Values are np and npq @ Finally, we can simulate independent random values from this distribution using \texttt{rbinom}. To simulate $N=50$ random, Binomial random variables: <>= N<-50 Ysamp<-rbinom(N,n,p) Ysamp @ \item \textbf{The Geometric distribution.} The geometric pmf $p(y)$ takes the form \[ p(y) = (1-p)^{y-1} p \qquad y = 1,2,\ldots, \] with $p(y) = 0$ for other values of $y$. This pmf can be computed using the \texttt{dgeom} function in \texttt{R}, however, \texttt{R} uses the alternative form of the pmf \[ p(y) = (1-p)^{y} p \qquad y = 0,1,2,\ldots, \] so we need to subtract one from the argument of \texttt{dgeom} <>= p<-0.2;y<-1:20 #Set the parameters p.y<-dgeom(y-1,p) #Compute the pmf par(mar=c(4,4,1,2)) #Set up the plot plot(y,p.y,pch=19,ylab='p(y)',cex=0.75) #Produce the plot print(cbind(y,p.y)) #Print the numerical values of y and p(y) @ We can use the function \texttt{pgeom} to compute the cdf $F(y)$: for $y = 1,2,3,\ldots$, \[ F(y) = P(Y \leq y) = \sum_{j=1}^y p(j) = \sum_{j=1}^{y} (1-p)^{y-1} p = 1 - (1-p)^y \] <>= F.y<-pgeom(y-1,p) #Compute the cdf par(mar=c(4,4,1,2)) plot(stepfun(y,c(0,F.y)),xlab='y',ylab='F(y)',pch=19,cex=0.5,vertical=FALSE,main='') @ We can compute the expectation and variance by using the \texttt{sum} function. Here, we can only approximate the infinite sum by taking the first 10000 (say) terms <>= y<-1:10000 #Take the first 10000 y values p.y<-dgeom(y-1,p) Expectation<-sum(y*p.y) #Compute the expectation Variance<-sum((y-Expectation)^2*p.y) #Compute the variance c(Expectation,Variance) #Values are 1/p and (1-p)/p^2 @ Finally, we can simulate independent random values from this distribution using \texttt{rgeom}. To simulate $N=50$ random, Geometric random variables: we must remember to add one to the generated values due to the alternative pmf that $\texttt{R}$ uses. <>= N<-50 Ysamp<-rgeom(N,p)+1 Ysamp @ \item \textbf{The Negative Binomial distribution.} The negative binomial pmf $p(y)$ takes the form \[ p(y) = \dbinom{y-1}{r-1} p^r (1-p)^{y-r} \qquad y = r,r+1,r+2,\ldots, \] with $p(y) = 0$ for other values of $y$, for some fixed positive integer $r$. This pmf can be computed using the \texttt{dnbinom} function in \texttt{R}, however, \texttt{R} uses the alternative form of the pmf \[ p(y) = \dbinom{y+r-1}{r-1} p^r (1-p)^{y} \qquad y = 0,1,2,\ldots, \] so we need to subtract $r$ from the argument of \texttt{dnbinom}. Below, we use $r=3$; <>= p<-0.2;r<-3;y<-r+c(0:20) #Set the parameters p.y<-dnbinom(y-1,r,p) #Compute the pmf par(mar=c(4,4,1,2)) #Set up the plot plot(y,p.y,pch=19,ylab='p(y)',cex=0.75,xlim=range(0,max(y)),ylim=range(0,0.2)) print(cbind(y,p.y)) #Print the numerical values of y and p(y) @ We can use the function \texttt{pnbinom} to compute the cdf $F(y)$: <>= F.y<-pnbinom(y-r,r,p) #Compute the cdf par(mar=c(4,4,1,2)) plot(stepfun(y,c(0,F.y)),xlab='y',ylab='F(y)',xlim=range(0,max(y)),ylim=range(0,1), pch=19,cex=0.5,vertical=FALSE,main='') @ We can compute the expectation and variance by using the \texttt{sum} function. Here, we can only approximate the infinite sum by taking the first 10000 (say) terms <>= y<-r+c(0:10000) #Take the first 10000 y values p.y<-dnbinom(y-r,r,p) Expectation<-sum(y*p.y) #Compute the expectation Variance<-sum((y-Expectation)^2*p.y) #Compute the variance c(Expectation,Variance) #Values are r/p and r(1-p)/p^2 @ Finally, we can simulate independent random values from this distribution using \texttt{rnbinom}. To simulate $N=50$ random, negative binomial random variables: we must remember to add $r$ to the generated values due to the alternative pmf that $\texttt{R}$ uses. <>= N<-50;Ysamp<-rnbinom(N,r,p)+r; Ysamp @ The experimental context that defines the use of the negative binomial distribution is that the random variable records the number of independent binary trials required to obtain the $r$th `success'. Therefore, we can regard $Y$ as being the sum of independent geometric random variables \[ Y = \sum_{i=1}^r Y_i \] where $Y_i$ is the number of trials between the $(i-1)$st and the $i$th success. For example, we can regard the sequence \[ 001000001001 \] by breaking it down as the three sub-sequences \[ 001 \quad + \quad 000001 \quad + \quad 001. \] We verify this in the following simulation of $N=100000$ negative binomial random variables with $r=3$ and $p=0.3$: <>= N<-100000 Yi.mat<-matrix(rgeom(N*3,p),ncol=3) #Simulate an N x 3 matrix of Geometric(0.2) r.v.s Y<-rowSums(Yi.mat)+3 #Form the vector of row sums of the matrix. Sim.tab<-cbind(1:25,dnbinom(1:25-r,r,p),tabulate(Y,nbins=25)/N) colnames(Sim.tab)<-c('y','p(y)','Rel. Freq.') Sim.tab @ In this table, the column \texttt{Rel. Freq.} is the relative frequency of times (out of 1000000 simulated values) that that value of $y$ was observed. We could expect roughly a proportion $p(y)$ of simulated values to take the value $y$: that is, the second and third columns should match. \end{enumerate} \end{document}