library(ggplot2) # Generating samples from different distributions rnorm(10) x <- rnorm(1000) qplot(x, geom="histogram", fill=I("black")) x <- rnorm(1000, mean=10, sd=5) qplot(x, geom="histogram", fill=I("black")) rpois(10, 10) x <- rt(500, 2) qplot(x, geom="histogram", fill=I("black")) x <- rcauchy(100, 2) qplot(x, geom="histogram", fill=I("black")) x <- rgamma(1000, 0.5) qplot(x, geom="histogram", fill=I("black")) qplot(x, geom="density", fill=I("black")) x <- rgamma(1000, 0.1, fill=I("black")) qplot(x, geom="density") x <- rgamma(1000, 1) qplot(x, geom="density", fill=I("black")) # Your turn qplot(rnorm(100), geom="histogram", fill=I("black"), binwidth=0.5) qplot(rnorm(100), geom="density", fill=I("black")) qplot(rnorm(50, 10, 5), geom="density", fill=I("black")) qplot(rpois(1000, 50), geom="histogram", fill=I("black"), binwidth=5) qplot(rbeta(100, 0.1, 0.1), geom="density", fill=I("black")) qplot(runif(30), geom="histogram", fill=I("black"), binwidth=0.1) # Repetition replicate(10, mean(rnorm(100))) qplot(replicate(100, mean(rnorm(10))),geom="histogram", fill=I("black"),binwidth=0.2, xlim=c(-1,1)) qplot(replicate(100, mean(rnorm(100))),geom="histogram", fill=I("black"), binwidth=0.05, xlim=c(-1,1)) # Your turn qplot(replicate(1000, mean(runif(2, 0, 10))),geom="histogram", fill=I("black"), binwidth=0.5, xlim=c(0,10)) qplot(replicate(1000, mean(runif(10, 0, 10))),geom="histogram", fill=I("black"), binwidth=0.5, xlim=c(0,10)) qplot(replicate(1000, mean(runif(100, 0, 10))),geom="histogram", fill=I("black"), binwidth=0.2, xlim=c(0,10)) # Central Limit Theorem # Functions clt <- function(n, m) { x <- replicate(n, mean(runif(m))) qplot(x, geom="histogram", fill=I("black")) } clt(1000, 2) clt(1000, 10) clt(1000, 100) # Adding flexibility clt <- function(n, m, r) { x <- replicate(n, mean(r(m))) qplot(x, geom="histogram", fill=I("black")) } clt(1000, 2, runif) clt(1000, 2, rexp) clt(1000, 20, rexp) clt(1000, 1, rpois) qplot(replicate(1000, mean(rpois(2, 1))), geom="histogram", fill=I("black")) clt(1000, 2, function(m) rpois(m, lambda=1)) clt <- function(n, m, r, ...) { f <- function() mean(r(m, ...)) x <- replicate(n, f()) qplot(x, geom="histogram", fill=I("black"), ...) } clt(n=1000, m=2, r=rnorm, mean = 10, sd = 10) clt(n=1000, m=2, r=rpois, lambda=1) clt(n=1000, m=2, r=rpois, lambda=10) clt(n=1000, m=2, r=rpois, lambda=0.5) # Exploring t-tests ?t.test t.test(rnorm(20), rnorm(20)) str(t.test(rnorm(20), rnorm(20))) t.test(rnorm(20), rnorm(20))\$p.value tt.check <- function(n) { t.test(rnorm(n), rnorm(n))\$p.value } qplot(replicate(5000, tt.check(10)), geom="histogram", fill=I("black"), binwidth=0.1) qplot(replicate(5000, tt.check(100)), geom="histogram", fill=I("black"), binwidth=0.1) qplot(replicate(5000, tt.check(1000)), geom="histogram", fill=I("black"), binwidth=0.1) tt.check <- function(n, mn=1) { t.test(rnorm(n), rnorm(n, mn))\$p.value } qplot(replicate(5000, tt.check(10)), geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) qplot(replicate(5000, tt.check(20)), geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) qplot(replicate(5000, tt.check(30)), geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) qplot(replicate(5000, tt.check(10, 0.5)), geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) qplot(replicate(5000, tt.check(20, 0.5)), geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) qplot(replicate(5000, tt.check(30, 0.5)), geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) # Check proportion of times: alpha=0.05 # Both means are equal, but reject the null hypothesis x<-replicate(5000, tt.check(10, 0)) length(x[x<0.05])/5000 x<-replicate(5000, tt.check(20, 0)) length(x[x<0.05])/5000 x<-replicate(5000, tt.check(30, 0)) length(x[x<0.05])/5000 # Both means are different, but fail to reject the null hypothesis x<-replicate(5000, tt.check(10)) length(x[x>0.05])/5000 x<-replicate(5000, tt.check(20)) length(x[x>0.05])/5000 x<-replicate(5000, tt.check(30)) length(x[x>0.05])/5000 x<-replicate(5000, tt.check(10, 0.5)) length(x[x>0.05])/5000 x<-replicate(5000, tt.check(20, 0.5)) length(x[x>0.05])/5000 x<-replicate(5000, tt.check(30, 0.5)) length(x[x>0.05])/5000 # What happens with samples from different different distributions tt.check.dist <- function(n, r1, r2) { t.test(r1(n), r2(n))\$p.value } # Uniform x<-replicate(5000,tt.check.dist(10, r1=runif, r2=runif)) length(x[x<0.05])/5000 x<-replicate(5000,tt.check.dist(20, r1=runif, r2=runif)) length(x[x<0.05])/5000 x<-replicate(5000,tt.check.dist(30, r1=runif, r2=runif)) length(x[x<0.05])/5000 qplot(x, geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) # Exponential x<-replicate(5000,tt.check.dist(10, r1=rexp, r2=rexp)) length(x[x<0.05])/5000 x<-replicate(5000,tt.check.dist(20, r1=rexp, r2=rexp)) length(x[x<0.05])/5000 x<-replicate(5000,tt.check.dist(30, r1=rexp, r2=rexp)) length(x[x<0.05])/5000 qplot(x, geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) # t, df=2, df=30 x<-replicate(5000, tt.check.dist(10, r1=function(n) rt(n,2), r2=function(n) rt(n,2))) length(x[x<0.05])/5000 x<-replicate(5000, tt.check.dist(20, r1=function(n) rt(n,2), r2=function(n) rt(n,2))) length(x[x<0.05])/5000 x<-replicate(5000, tt.check.dist(30, r1=function(n) rt(n,2), r2=function(n) rt(n,2))) length(x[x<0.05])/5000 qplot(x, geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) x<-replicate(5000, tt.check.dist(10, r1=function(n) rt(n,30), r2=function(n) rt(n,30))) length(x[x<0.05])/5000 x<-replicate(5000, tt.check.dist(20, r1=function(n) rt(n,30), r2=function(n) rt(n,30))) length(x[x<0.05])/5000 x<-replicate(5000, tt.check.dist(30, r1=function(n) rt(n,30), r2=function(n) rt(n,30))) length(x[x<0.05])/5000 qplot(x, geom="histogram", fill=I("black"), binwidth=0.1, xlim=c(0,1)) # Bootstrap ?galaxies summary(galaxies) qplot(galaxies, geom="density", fill=I("black")) median(sample(galaxies, replace=T)) x<-replicate(5000, median(sample(galaxies, replace=T))) mean(x) sd(x) qplot(x, geom="density", fill=I("black"), xlim=range(galaxies)) qplot(x, geom="histogram", fill=I("black")) # Your turn qplot(diamonds\$price, geom="density", fill=I("black")) x<-replicate(1000, var(sample(diamonds\$price, replace=T))) mean(x) sd(x) qplot(x, geom="density", fill=I("black")) x<-replicate(1000, var(sample(diamonds\$price, 10, replace=T))) mean(x) sd(x) qplot(x, geom="density", fill=I("black")) # Permutation x<-cbind(rnorm(10), rnorm(10)) t.test(x[,1],x[,2])\$statistic y<-as.vector(x) lb<-rep(c(1,2), rep(nrow(x),2)) y<-cbind(y,sample(lb)) colnames(y)<-c("x","label") t.test(y[lb==1],y[lb==2])\$p.value qplot(x,data=data.frame(y), geom="histogram", facets=label~.) tt.check.permute <- function(x) { y<-as.vector(x) lb<-rep(c(1,2), rep(nrow(x),2)) y<-cbind(y,sample(lb)) t.test(y[lb==1],y[lb==2])\$statistic } x<-cbind(rnorm(10), rnorm(10)) y<-as.vector(x) lb<-rep(c(1,2), rep(nrow(x),2)) y<-cbind(y,lb) colnames(y)<-c("x","label") qplot(x,data=data.frame(y), geom="histogram", fill=I("black"), facets=label~.) z<-replicate(1000, tt.check.permute(x)) qplot(z, geom="histogram", fill=I("black")) t<-abs(t.test(x[,1],x[,2])\$statistic) length(z[z<(-t) | z>t])/1000 x<-cbind(rnorm(100), rnorm(100)) y<-as.vector(x) lb<-rep(c(1,2), rep(nrow(x),2)) y<-cbind(y,lb) colnames(y)<-c("x","label") qplot(x,data=data.frame(y), geom="histogram", fill=I("black"), facets=label~.) z<-replicate(1000, tt.check.permute(x)) qplot(z, geom="histogram", fill=I("black")) t<-abs(t.test(x[,1],x[,2])\$statistic) length(z[z<(-t) | z>t])/1000 x<-cbind(rnorm(100), rnorm(100, 0.5)) y<-as.vector(x) lb<-rep(c(1,2), rep(nrow(x),2)) y<-cbind(y,lb) colnames(y)<-c("x","label") qplot(x,data=data.frame(y), geom="histogram", fill=I("black"), facets=label~.) z<-replicate(1000, tt.check.permute(x)) qplot(z, geom="histogram", fill=I("black")) t<-abs(t.test(x[,1],x[,2])\$statistic) length(z[z<(-t) | z>t])/1000 # Graphically qplot(carat, price, data=diamonds, geom=c("point","smooth"), method="lm") for (i in 1:5) { p <- qplot(carat, sample(price), data=diamonds, geom=c("point","smooth"), method="lm") print(p) readline("Ready to continue?") }