# Chihara and Hesterberg, chap 4 ### Chapter 4: Sampling Distributions #--------------------------------------------- # Example 4.1: Sampling distribution from Bin(size=10, p=1/2) # toss a fair coin 10 times and note the proportion of heads, pHat nHeads <- rbinom(n=1, size=10, p=1/2) pHat <- nHeads/10 # one experiment : record the proportion of heads in 10 coin flips # repeat the experiment 1000 times my.means<-numeric(1000) #set.seed(300) for (i in 1:1000) { x <- rbinom(n=1, size=10, p=1/2) my.means[i]<- mean(x)/10 } table(my.means) plot(table(my.means), col="orange3", xlab="pHat", ylab="Count", main="Distribution of pHat after 1000 sets of 10 tosses") hist(my.means, main="Simulated sampling distribution", xlim=c(0,1), ylim=c(0,2.8), col="palegoldenrod", xlab="sample means", prob=TRUE) curve(dnorm(x, mean=1/2, sd=sqrt(1/40)), from=0, to=1, col="indianred2", add=TRUE) dev.new() qqnorm(my.means, col="lightskyblue2") qqline(my.means, col="blue") # mean == 1/2 mean(my.means) # sd = sqrt(p*q/n) sd(my.means) #--------------------------------------------- # Example 4.2: Sampling distribution from Exp(1/15) my.means<-numeric(1000) #set.seed(300) for (i in 1:1000) { x<-rexp(100, 1/15) my.means[i]<- mean(x) } hist(my.means, main="Simulated sampling distribution from Exp(1/15)", col="steelblue2", xlab="means") dev.new() qqnorm(my.means, col="yellowgreen") qqline(my.means, col="violetred4") mean(my.means) sd(my.means) #---------------------------------------------------- # Example 4.3: Sampling Distribution from Unif[0,1] my.max<-numeric(1000) #set.seed(100) for (i in 1:1000) { y <- runif(12) #draw random sample of size 12 my.max[i] <- max(y) #find max, save in i_th position of my.max } hist(my.max, col="seagreen3", main="Sampling Distribution of Maxima from Unif[0,1]", xlab="maximums") # To create a histogram with a density curve imposed # scale bars to have area one with prob=TRUE option hist(my.max, col="seashell2", main="Sampling Distribution of Maxima from Unif[0,1]", xlab = "maximums", prob=TRUE) #add pdf to histogram curve(12*x^{11}, col="blue", add=TRUE) #------------------------------------------------ # Example 4.6 # Sampling distribution simulation # Sample of size 30 from gamma r=5, lambda=2 #set.seed(10) my.means<-numeric(1000) for (i in 1:1000) { x <- rgamma(30, shape=5, rate=2) my.means[i] <- mean(x) } hist(my.means, col="turquoise2", main="Distribution of means from Gamma(r=5, lambda=2)") dev.new() qqnorm(my.means, col="wheat4") qqline(my.means, col="yellow2") mean(my.means) sd(my.means) sum(my.means > 3)/1000 #alternatively mean(my.means > 3) #----------------------------------------------