﻿ ## ## Mean ## mu <- 178 ## ## Standard deviation ## sigma <- 12 ## ## Sample size ## n <- 10 ## ## Simulate normally distributed X_i ## x <- rnorm(n=n, mean=mu, sd=sigma) ## ## See the realizations ## x ## ## Empirical density ## hist(x, prob=TRUE, col='blue') ## ## Find the sample mean ## mean(x) ## ## Find the sample variance ## var(x) ## ## Repeat the simulated sampling many times ## mat <- replicate(100, rnorm(n=n, mean=mu, sd=sigma)) ## ## Find the sample mean for each of them ## xbar <- apply(mat, 2, mean) ## ## Now we have many realizations of the sample mean ## xbar ## ## See their distribution ## hist(xbar, prob=TRUE, col='blue') ## ## There mean ## mean(xbar) ## ## and sample variance ## var(xbar) x=seq(-4,4, by=0.01) plot(x,dt(x,9), type="l", col=2) lines(x,dnorm(x), type="l") text(2.5,0.3,"Black: standard normal") text(3,0.1,"Red: t(9)", col = 2) x=seq(-4,4, by=0.01) plot(x,dt(x,9), type="l", col=2) lines(x,dnorm(x), type="l") text(2.5,0.3,"Black: standard normal") text(3,0.1,"Red: t(9)", col = 2) polygon(c(2,seq(2,4, by=0.01),4,2),c(0,dt(seq(2,4, by=0.01),9),0, 0),col="pink") text(3.5,0.04, "P(T>2)=0.038",col=2) x=seq(-4,4, by=0.01) plot(x,dt(x,9), type="l", col=2) lines(x,dnorm(x), type="l") text(2.5,0.3,"Black: standard normal") text(2.8,0.1,"Red: t(9)", col = 2) polygon(c(2,seq(2,4, by=0.01),4,2),c(0,dnorm(seq(2,4, by=0.01)),0, 0),col="grey") text(3.3,0.04, "P(Z>2)=0.023") ## The t-quantiles for n=10: qt(0.975,9) qt(0.995,9) x <- c(168,161,167,179,184,166,198,187,191,179) t.test(x,conf.level=0.99) par(mar = c(4.5, 4.5, 3.5, 0.5)) n=1 k=1000 u=matrix(runif(k*n),ncol=n) hist(apply(u,1,mean),col="blue",main="n=1",xlab="Means") par(mar = c(4.5, 4.5, 3.5, 0.5)) n=2 k=1000 u=matrix(runif(k*n),ncol=n) hist(apply(u,1,mean),col="blue",main="n=2",xlab="Means") par(mar = c(4.5, 4.5, 3.5, 0.5)) n=6 k=1000 u=matrix(runif(k*n),ncol=n) hist(apply(u,1,mean),col="blue",main="n=6",xlab="Means") par(mar = c(4.5, 4.5, 3.5, 0.5)) n=30 k=1000 u=matrix(runif(k*n),ncol=n) hist(apply(u,1,mean),col="blue",main="n=30",xlab="Means", nclass=15) x=seq(-3.5,3.5,by=0.01) plot(x,dt(x,9),type="l",xaxt="n",xlab="",ylab="",yaxt="n") polygon(c(qt(0.975,9),seq(qt(0.975,9),4, by=0.01),4, qt(0.975,9)),c(0,dt(seq(qt(0.975,9),4, by=0.01),9),0, 0),col="pink") polygon(c(-qt(0.975,9),seq(-qt(0.975,9),-4, by=-0.01),-4,-qt(0.975,9)),c(0,dt(seq(-qt(0.975,9),-4, by=-0.01),9),0,0),col="pink") text(0,0.1,expression(P(group("|",frac(bar(X)-mu,S/sqrt(n)),"|")