\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{Sampling Distributions: Examples} \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 A \textit{sampling distribution} is simply the probability distribution of a summary statistic random variable: we consider \begin{itemize} \item a \textit{random sample} of (independent and identically distributed) random variables $X_1,\ldots,X_n$ having some pmf $p(x)$ or pdf $f(x)$ \item a summary statistic random variable, for example $\Xbar$, derived from $X_1,\ldots,X_n$ \item the probability distribution of the summary statistic, calculated from the pmf or pdf of the original variables. \end{itemize} Note: these mathematical calculations are performed \textbf{without} considering any observed data (which would be the numerical values $x_1,\ldots,x_n$). The calculation is allowing us to specify how the summary statistic should behave if we make assumptions about the distributions of the original random variables. \bigskip <>= x<-2:5 px<-c(0.05,0.15,0.6,0.2) mu<-sum(x*px) sigmasq<-sum(x^2*px)-sum(x*px)^2 sigma<-sqrt(sigmasq) @ \textbf{Example:} Suppose random variable $X$ records the number of courses in which a typical student at McGill is registered in a given term, and suppose that the probability distribution of $X$ is specified by the following pmf \begin{table}[ht] \centering \begin{tabular}{|c|c|c|c|c|} \hline % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ... $x$ & 2 & 3 & 4 & 5 \\ \hline $p(x)$ & 0.05 & 0.15 & 0.6 & 0.2 \\ \hline \end{tabular} \end{table} \noindent That is, the probability that a student is registered in four courses is 0.6 etc. We can compute that the expected number of courses in which a student is registered is \[ E(X) = \sum_{x=1}^4 x p(x) = (\Sexpr{x[1]} \times \Sexpr{px[1]}) + (\Sexpr{x[2]} \times \Sexpr{px[2]}) + (\Sexpr{x[3]} \times \Sexpr{px[3]}) + (\Sexpr{x[4]} \times \Sexpr{px[4]}) = \Sexpr{sum(x*px)} \] and for the variance \[ E(X^2) = \sum_{x=1}^4 x^2 p(x) = (\Sexpr{x[1]}^2 \times \Sexpr{px[1]}) + (\Sexpr{x[2]}^2 \times \Sexpr{px[2]}) + (\Sexpr{x[3]}^2 \times \Sexpr{px[3]}) + (\Sexpr{x[4]}^2 \times \Sexpr{px[4]}) = \Sexpr{sum(x^2*px)} \] so \[ Var(X) = E(X^2) - \{E(X)\}^2 = \Sexpr{sum(x^2*px)} - \Sexpr{sum(x*px)^2} = \Sexpr{sum(x^2*px)-sum(x*px)^2} \] so the standard deviation is $\sqrt{\Sexpr{sigmasq}} = \Sexpr{myrnd(sigma,4)}$. \bigskip Suppose we intend to consider $n=3$ students, and will record the number of courses in which they are registered as random variables $X_1$, $X_2$ and $X_3$. Before we record the number of courses for each, we can compute the probability distribution of the sample mean random variable \[ \Xbar = \frac{X_1+X_2+X_3}{3} \] by considering all possible values of the three variables. There are 64 possible configurations, and these are listed in the following table: <>= x1<-rep(2:5,each=16) x2<-rep(rep(2:5,each=4),2) x3<-rep(2:5,16) xposs<-cbind(x1,x2,x3) xsum<-rowSums(xposs) xbar<-xsum/3 xprob<-apply(xposs,1,function(x){return(prod(px[x-1]))}) xconf<-cbind(xposs,xsum,xbar,xprob) colnames(xconf)<-c('1','2','3','Sum','Mean','Prob.') library(kableExtra) xptab<-cbind(xconf[1:32,],xconf[33:64,]) kbl(xptab, "latex",longtable = T, booktabs = T, caption = "Possible configurations") %>% add_header_above(c("Configurations 1 to 32" = 6, "Configurations 33 to 64" = 6)) %>% kable_styling(latex_options = c("repeat_header")) @ The column `Prob' is the probability of obtaining that configuration, which we compute using independence: \[ P(X_1 = x_1, X_2 = x_2, X_3 = x_3) = P(X_1 = x_1) P(X_2 = x_2) P(X_3 = x_3) \] for all possible values of $x_1$, $x_2$, $x_3$. From this table, we deduce that the sample mean random variable can take $\Sexpr{length(table(xbar))}$ possible values, namely \[ \calX = \{\Sexpr{myrnd(as.numeric(names(table(xbar))),4)}\} \] and we can compute the probability of each possibility simply by adding the probabilities across configurations that give the same mean value. <>= df<-data.frame(xconf) suppressMessages(library(tidyverse)) df.sum<-df %>% group_by(Mean) %>% summarize(totalprob=sum(Prob.)) df.sum<-as.data.frame(df.sum) names(df.sum)[2]<-'Probability' kbl(df.sum, "latex",longtable = T, booktabs = T, caption = "Distribution of Mean random variable") %>% kable_styling(latex_options = c("repeat_header")) @ That is, we have that the pmf of $\Xbar$ is given by the probabilities in this table. We can denote this pmf $p_{\Xbar}(x)$ (with the subscript indicating which random variable we are dealing with), and plot it in the usual way. <>= par(mar=c(4,4,2,0)) plot(df.sum[,1],df.sum[,2],pch=19,ylim=range(0,0.5), xaxt='n',xlab='x',ylab=expression(paste('PMF of ',bar(X)))) axis(1,at=myrnd(df.sum[,1],3)) @ \noindent We can treat this as a standard pmf on the set of values $\calX$: we can compute the expectation \[ E(\Xbar) = \sum_{x \in \calX} x p_{\Xbar}(x) = \Sexpr{sum(df.sum[,1]*df.sum[,2])} \] which is \textbf{identical} to the original expectation $E(X)$. For the variance, however, we have by an identical method of calculation that \[ Var(\Xbar) = \Sexpr{sum((df.sum[,1]-sum(df.sum[,1]*df.sum[,2]))^2*df.sum[,2])} = \frac{Var(X)}{3}. \] What this result tells is that if we assume the original discrete distribution, and we assume that our three students are drawn from this population, then the sample mean that we will get by averaging the number of courses they take will \begin{itemize} \item take values in $\calX$, \item have probability distribution described by pmf $p_{\Xbar}$ \item have expectation equal to the expected number from the original distribution \item have variance equal to the variance from the original distribution divided by three. \end{itemize} We can test this in simulation: suppose we plan to obtain data on 10000 batches of three students. We can use the \texttt{sample} function: <>= set.seed(2) X.sample<-matrix(sample(x,px,size=30000,rep=T),nrow=10000) Xmean<-apply(X.sample,1,mean) df.sum$rfreq<-as.numeric(table(Xmean)/10000) df.sum mean(Xmean) #observed sample mean var(Xmean) #observed sample variance @ The relative frequency of obtaining each of the possible 10 values of the mean matches closely the computed pmf, and the sample mean and sample variance are close to their theoretical population versions. \bigskip Note: Here we have computed for samples of $n=3$. We can repeat the same calculation for any $n$; however, there will be $4^n$ possible configurations, and this number grows rapidly, so the exact computation method outlined above is difficult to implement. \pagebreak \textbf{Example:} Recall the body temperature data from earlier in the course; these are $n=136$ measurements of body temperature (in Celsius) <>= BodyTemp<-scan('https://www.math.mcgill.ca/dstephens/203/Data/BodyTemp.txt') n<-length(BodyTemp) c(mean(BodyTemp),var(BodyTemp),sd(BodyTemp)) @ For these observed data, we have the sample summaries \[ \xbar = \Sexpr{mean(BodyTemp)} \qquad s^2 = \Sexpr{var(BodyTemp)} \qquad s = \Sexpr{sd(BodyTemp)} \] Suppose we plan to carry out a similar measurement experiment by collecting data on $n=136$ patients in each of ten hospitals, where we believe that the body temperature of all patients will follow a $Normal(\mu, \sigma^2)$ distribution with $\mu = 37$ (Celsius) and $\sigma = 0.40$. We might observe the following sample means once we have collected the data: <>= X<-matrix(rnorm(n*10,37,0.40),nrow=10) Xmean<-apply(X,1,mean) Xmean @ As we see, there is some variation in the observed sample means across hospitals, which merely reflects sampling variability. We can prove mathematically that the sample mean random variable has a distribution which is also Normal, that is \[ \Xbar \sim Normal \left( \mu, \frac{\sigma^2}{n} \right) \] and we have thus characterized the probability distribution of $\Xbar$. Note that if we consider the $Z$-transform method, we have that \[ Z = \frac{\Xbar - \mu}{\sigma/\sqrt{n}} = \sqrt{n}\frac{(\Xbar-\mu)}{\sigma} \sim Normal(0,1) \] \end{document}