# Why poor Gosset went to graduate school: When all goes well when brewing Guiness, a sample of the beer will have X amount of "ick" in it, where X= N(mu,sigma) for some known mu and unknown sigma. We can view this as our Null hypothesis. When thing go wrong mu1>mu. If mu1 tends to be a lot bigger than mu when things go wrong then even a very small number of samples might reveal this. If we new sigma even a single sample might do. Since sigma was mysterious, Gosset used 4 sample and estimated the true value of the mean via (X1+X2+X3+X4)/4. Under the Null hypothesis is N(mu0,sigma/sqrt(4)). Imagine, he wanted a test such that when things went right he only claimed they went wrong 5 percent of the time. So he concluded that he should use the rule that if muhat-mu> (1.645)*SE = (1.645)*s/sqrt(4) then he would reject the Null. Let's run it.... ErrorRate=0.05 mu=15 sigma=3 n=4 X=rnorm(n,mu,sigma) Xbar=mean(X) s=sd(X) Xbar s zstar=abs(qnorm(ErrorRate,0,1)) SE=s/sqrt(n) Xbar-zstar*SE ################## #Lets do it many times PreGosset<-function(Trials,n,mu,sig,ErrorRate){ H=c() for (i in 1:Trials){ X=rnorm(n,mu,sigma) Xbar=mean(X) s=sd(X) SE=s/sqrt(n) zstar=abs(qnorm(ErrorRate,0,1)) H<-rbind(H,c(Xbar-zstar*SE,Xbar,Xbar+zstar*SE)) } H } Trials=400 n=4 mu=15 sigma=3 ErrorRate=0.05 Confs=PreGosset(Trials,n,mu,sigma,ErrorRate) left=min(Confs[,1])-.3 right=max(Confs[,3])+.3 Count=0 plot(c(left,left,right,right),c(Trials+1,0,0,Trials+1),type="l",main=paste(Trials,"Pre-Gosset Trial with error rate",ErrorRate), ylab="Trials",xlab="ick") for (j in 1:Trials){ points(Confs[j,c(1,3)],c(j,j),type="l",col="red") #points(Confs[j,2],c(j),pch=18) points(Confs[j,1],c(j),pch=18) {if (Confs[j,1]>mu) Count=Count+1} #points(Confs[j,1],c(j+1),pch=40) #points(Confs[j,3],c(j+1),pch=41) } Count points(c(mu,mu),c(0,Trials+1),type="l",col="blue") # After he claims a sample is too icky, the actually verifies his claim. Why was Gosset nearly fired? In other words, what fraction the non-icky beer batches did Gosset have the lab test? What happened?!? # Answer: Using s to approximate sigma created the problem. Namely, if the mean is high (or low) then the samples are likely to all be high and hence their variation small. This can be compensated for by using a larger "zstar" than you'd expect, and it's called a "tstar". It depends on the number of samples via what is called the degrees of freedom, denoted df. df=n-1 tstar=abs(qt(ErrorRate,df)) tstar ################## # Before explaining tstar further, let us verify that it indeed works: Gosset<-function(Trials,n,mu,sig,ErrorRate){ H=c() for (i in 1:Trials){ X=rnorm(n,mu,sigma) Xbar=mean(X) s=sd(X) SE=s/sqrt(n) df=n-1 tstar=abs(qt(ErrorRate,df)) H<-rbind(H,c(Xbar-tstar*SE,Xbar,Xbar+tstar*SE)) } H } Trials=400 n=4 mu=15 sigma=3 ErrorRate=0.05 Confs=Gosset(Trials,n,mu,sigma,ErrorRate) left=min(Confs[,1])-.3 right=max(Confs[,3])+.3 Count=0 plot(c(left,left,right,right),c(Trials+1,0,0,Trials+1),type="l",main=paste(Trials,"Post-Gosset Trial with ErrorRate rate",ErrorRate), ylab="Trials",xlab="ick") for (j in 1:Trials){ points(Confs[j,c(1,3)],c(j,j),type="l",col="red") #points(Confs[j,2],c(j),pch=18) points(Confs[j,1],c(j),pch=18) {if (Confs[j,1]>mu) Count=Count+1} #points(Confs[j,1],c(j+1),pch=40) #points(Confs[j,3],c(j+1),pch=41) } Count points(c(mu,mu),c(0,Trials+1),type="l",col="blue") ################## # To understand this better we can turn this confidence interval view on its head and consider the hypothesis test. We ask "Under the Null, what is the distribution of t=(Xbar-mu0)/(s/sqrt(n))?" Well, this distribution is not normal and it has a funny name: the student-t distribution with df degrees of freedom. Let us empirically explore this. mu=15 sigma=3 n=4 N=10000 M=n*N X=rnorm(M,mu,sigma) dim(X)<-c(N,n) Xbar=apply(X,1,mean) s=apply(X,1,sd) T=(Xbar-mu)/(s/sqrt(n)) #ST=sort(T) CutOff=8 hist(T[abs(T)